# An Interior-Point Trust-Region Method for Nonsmooth Regularized Bound-
Constrained Optimization

Geoffroy Leconte GERAD and Department of Mathematics and Industrial
Engineering, Polytechnique Montréal, QC, Canada. E-mail:
<EMAIL_ADDRESS> Dominique Orban GERAD and Department of
Mathematics and Industrial Engineering, Polytechnique Montréal, QC, Canada.
E-mail<EMAIL_ADDRESS>Research partially supported by an NSERC
Discovery Grant.

###### Abstract

We develop an interior-point method for nonsmooth regularized bound-
constrained optimization problems. Our method consists of iteratively solving
a sequence of unconstrained nonsmooth barrier subproblems. We use a variant of
the proximal quasi-Newton trust-region algorithm TR of Aravkin et al. [2] to
solve the barrier subproblems, with additional assumptions inspired from well-
known smooth interior-point trust-region methods. We show global convergence
of our algorithm with respect to the criticality measure of Aravkin et al.
[2]. Under an additional assumption linked to the convexity of the nonsmooth
term in the objective, we present an alternative interior-point algorithm with
a slightly modified criticality measure, which performs better in practice.
Numerical experiments show that our algorithm performs better than the trust-
region method TR, the trust-region method with diagonal hessian approximations
TRDH of Leconte and Orban [17], and the quadratic regularization method R2 of
Aravkin et al. [2] for two out of four tested bound-constrained problems. On
those two problems, our algorithm obtains smaller objective values than the
other solvers using fewer objective and gradient evaluations. On the two other
problems, it performs similarly to TR, R2 and TRDH.

###### keywords:

Regularized optimization, nonsmooth optimization, nonconvex optimization,
bound constraints, proximal gradient method, barrier method.

49J52, 65K10, 90C53, 90C56,

##  1 Introduction

We consider the problem

(1)
$\mathop{\textup{minimize}}_{x\in\mathds{R}^{n}}f(x)+h(x)\quad\mathop{\textup{subject
to}}\ x\geq 0,$

where $f:\mathds{R}^{n}\to\mathds{R}$ has Lipschitz-continuous gradient with
constant $L_{f}\geq 0$ and $h:\mathds{R}^{n}\to\mathds{R}\cup\\{+\infty\\}$ is
proper and lower semi-continuous. Both $f$ and $h$ may be nonconvex. $h$ is
often considered as a regularization function used to favor solutions with
desirable properties, such as sparsity.

Problems such as (1) are classically solved with a variant of the proximal
gradient method [19]. The proximal quasi-Newton trust-region algorithm of
Aravkin et al. [2], referred to as TR, can be extended to box constraints,
provided that the proximal operator of $h+\chi(\cdot\mid[\ell,u])$, where
$\ell<u$ componentwise and $\chi$ is the indicator of the box $[\ell,u]$, can
be computed efficiently. Leconte and Orban [17] present a variant of TR named
TRDH that also supports box constraints, and uses diagonal quasi-Newton
approximations. For specific separable regularizers $h$, they provide a
closed-form solution of the trust-region subproblems with box constraints,
giving rise to what they defined as an indefinite proximal operator; a scaled
generalization of a proximal operator.

At each iteration $k$, we solve the barrier subproblem

(2) $\mathop{\textup{minimize}}_{x\in\mathds{R}^{n}}f(x)+\phi_{k}(x)+h(x),$

where $\phi_{k}$ is the logarithmic barrier function

(3)
$\phi_{k}(x)=\sum_{i=1}^{n}\phi_{k,i}(x),\qquad\phi_{k,i}(x):=-\mu_{k}\log(x_{i}),\
i=1,\ldots,n,$

and $\\{\mu_{k}\\}\searrow 0$. Each subproblem (2) is an unconstrained problem
if we consider that $-\log(x)=+\infty$ when $x\leq 0$. In Section 5.2, we
explain that under reasonable assumptions, we can solve solve (2) with a
modified version of Aravkin et al.’s TR algorithm, and we expect that the
solutions of (3) converge to a solution of (1) as $\mu_{k}\to 0$.

Our approach is sometimes referred to as a trust-region interior-point method,
or trust-region method for barrier functions. We refer the reader to [9,
Chapter $13$] for more information on the case where $h=0$. Our algorithm,
named RIPM (_Regularized Interior Proximal Method_), can be seen as a
generalization of those methods to solve (1).

An inconvenient of solving (2), induced by the logarithmic barrier function
(3), is that the smooth part of the subproblem $f+\phi_{k}$ does not have a
Lipschitz gradient, thus compromising the convergence properties of TR
established by Aravkin et al. [2]. Nevertheless, in our analysis, we establish
the convergence of the barrier subproblems using the update rules of Conn et
al. [9, Chapter $13.6.3$] for our trust-region model. We also show global
convergence of RIPM towards a first-order stationary point of (1) if the
trust-region radii and the step lengths used in proximal operator evaluations
are bounded away from zero, and the iterates generated by the algorithm remain
bounded. In Section 5.4, under a convexity assumption on the nonsmooth term,
we provide an alternative implementation of the outer iterations of RIPM where
we change the stopping criteria to improve numerical performance.

In addition, we implement a variant of RIPM named RIPMDH (_Regularized
Interior Proximal Method with Diagonal Hessian approximations_) that uses TRDH
to solve the barrier subproblems. We compare the performance of RIPM and
RIPMDH with TR, TRDH and R2, all available from RegularizedOptimization.jl
[5], on four bound-constrained problems. The first two problems are a
regularized box-constrained quadratic problem, and a sparse nonnegative matrix
factorization problem. These two problems require many TR and R2 iterations to
converge. RIPM and RIPMDH obtain smaller objective values than the other
solvers using fewer objective and gradient evaluations, which suggests that
they may be best suited to solve difficult bound-constrained nonsmooth
problems. The third problem is an inverse problem for finding the parameters
of a differential equation. RIPM and RIPMDH perform more objective and
gradient evaluations than TR, but RIPMDH performs the least amount of proximal
operator evaluations. The last problem is a regularized box-constrained basis-
pursuite denoise problem. RIPMDH exhibits similar performance to those of TR,
TRDH and R2 using a modification of some of its parameters.

#### Related research

Attouch and Wets [4] use $f=0$ and the nonsmooth barrier
$-\mu_{k}\sum_{i=1}^{n}\log(\min(\tfrac{1}{2},x_{i}))>0$. They use the theory
of epi-convergence to explain some convergence properties of barrier methods,
and in particular, that the objectives of the barrier subproblems epi-converge
to the objective of their initial constrained problem.

Chouzenoux et al. [7] use convex $f$ and $h$ with general inequality
constraints $c_{i}(x)\leq 0$ for $i\in\\{1,...,p\\}$, $p>0$ to solve large-
scale image processing problems. Their algorithm uses proximal gradient steps
to solve the barrier subproblem.

Bertocchi et al. [6] solve inverse problems $y=\mathcal{D}(H\bar{x})$, where
$y$ is some observed data, $\bar{x}$ the signal to determine, $H$ is a linear
observation operator, and $\mathcal{D}$ is an operator applying noise
perturbation. They solve the problem
$\mathop{\textup{minimize}}_{x\in\mathds{R}^{n}}f(Hx,y)+h(x)$ with constraints
$\ c_{i}(x)\geq 0,i\in\\{1,...,p\\}$, where $f$ is a convex function involving
the observations and the signal $x$, $f(\cdot;y)$ and $h$ are twice
differentiable, $h$ and $-c_{i}$ are convex (among other properties). They
compute the proximal operator of the barrier term, and use an interior-point
algorithm. They apply their algorithm to develop a neural network architecture
for image restoration.

De Marchi and Themelis [11] use the proximal gradient method to solve
nonsmooth regularized optimization problems such as (1) where $f$ has a
locally Lipschitz-continuous gradient, $h$ is continuous relative to its
domain and prox-bounded. In addition, the constraint $x\geq 0$ is replaced by
a more general constraint $c(x)\leq 0$, where $c$ has locally Lipschitz-
continuous Jacobian.

Shen et al. [23] present an active set proximal algorithm to solve (1) with
$h(x)=\lambda\|x\|_{1}$ for some $\lambda>0$, and where $\ell\leq x\leq u$
with $\ell<0<u$ instead of $x\geq 0$. They use a hybrid search direction based
upon a proximal gradient step for the active variables (i.e., the variables
that are at one the bounds of the constraints), and a Newton step for the
other variables.

#### Notation

For $v\in\mathds{R}^{n}$, $\|v\|$ denotes the Euclidean norm of $v$.
$\mathds{R}_{+}$ and $\mathds{R}_{++}$, denote, respectively, the sets of
positive and strictly positive real numbers, whereas $\mathds{R}_{+}^{n}$ and
$\mathds{R}_{++}^{n}$ denote the sets of vectors having all their components
in $\mathds{R}_{+}$ and $\mathds{R}_{++}$, respectively.

$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
denotes $\mathds{R}\cup\\{\pm\infty\\}$. The unit closed ball defined with the
$\ell_{\infty}$-norm and centered at the origin is $\mathds{B}$, and the ball
centered at the origin of radius $\Delta>0$ is $\Delta\mathds{B}$. If
$C\subseteq\mathds{R}^{n}$, the indicator of $C$ is $\chi(\cdot\mid
C):\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
defined by $\chi(x\mid C)=0$ if $x\in C$, and $\chi(x\mid C)=+\infty$ if
$x\not\in C$. For $y\in\mathds{R}^{n}$, the set $y+C$ is composed of all the
vectors $s\in\mathds{R}^{n}$ such that $s=y+x$ with $x\in C$.

Following the notation of Rockafellar and Wets [22], the set of all
subsequences of $\mathds{N}$ is denoted by $\mathcal{N}_{\infty}^{\\#}$, and
the set composed of the subsequences of $\mathds{N}$ containing all $k$ beyond
some $k_{0}$ is denoted by $\mathcal{N}_{\infty}$. For
$N\in\mathcal{N}_{\infty}^{\\#}$,
$\\{x_{k}\\}\underset{N}{\longrightarrow}\bar{x}$ indicates that the
subsequence $\\{x_{k}\\}_{k\in N}$ (which we may also write $\\{x_{k}\\}_{N}$
for conciseness) converges to $\bar{x}$.

$X$, $Z$ and $S$ (possibly with subscripts$~{}_{k}$,$~{}_{k,j}$
or$~{}_{k,j,1}$) denote the square diagonal matrices having $x$, $z$ and $s$
as diagonal elements, respectively.

##  2 Background

The following are standard variational analysis concepts—see, e.g., [22]. Let
$\phi:\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
and $\bar{x}\in\mathds{R}^{n}$ where $\phi$ is finite. The Fréchet
subdifferential of $\phi$ at $\bar{x}$ is the closed convex set
$\widehat{\partial}\phi(\bar{x})$ of elements $v\in\mathds{R}^{n}$ such that

$\liminf_{\begin{subarray}{c}x\to\bar{x}\\\
x\neq\bar{x}\end{subarray}}\frac{\phi(x)-\phi(\bar{x})-v^{T}(x-\bar{x})}{\|x-\bar{x}\|}\geq
0.$

The limiting subdifferential of $\phi$ at $\bar{x}$ is the closed, but not
necessarily convex, set $\partial\phi(\bar{x})$ of elements
$v\in\mathds{R}^{n}$ for which there exist $\\{x_{k}\\}\to\bar{x}$ and
$\\{v_{k}\\}\to v$ such that $\\{\phi(x_{k})\\}\to\phi(\bar{x})$ and
$v_{k}\in\widehat{\partial}\phi(x_{k})$ for all $k$. The inclusion
$\widehat{\partial}\phi(\bar{x})\subset\partial\phi(\bar{x})$ always holds.
Finally, the horizon subdifferential of $\phi$ at $\bar{x}$ is the closed, but
not necessarily convex, cone $\partial^{\infty}\phi(\bar{x})$ of elements
$v\in\mathds{R}^{n}$ for which there exist $\\{x_{k}\\}\to\bar{x}$,
$\\{v_{k}\\}$ and $\\{\lambda_{k}\\}\searrow 0$ such that
$\\{\phi(x_{k})\\}\to\phi(\bar{x})$, $v_{k}\in\widehat{\partial}\phi(x_{k})$
for all $k$, and $\\{\lambda_{k}v_{k}\\}\to v$.

If $\inf\phi>-\infty$, $\mathop{\textup{argmin}}\phi$ is the set of
$x\in\mathds{R}^{n}$ such that $\phi(x)=\inf\phi$. For $\epsilon>0$,
$\epsilon\text{-}\mathop{\textup{argmin}}\phi$ is the set of
$x\in\mathds{R}^{n}$ such that $\phi(x)\leq\inf\phi+\epsilon$.

If $C\subseteq\mathds{R}^{n}$ and $\bar{x}\in C$, the closed convex cone
$\widehat{N}_{C}(\bar{x}):=\widehat{\partial}\chi(\bar{x}\mid C)$ is the
regular normal cone to $C$ at $\bar{x}$. The closed cone
$N_{C}(\bar{x}):=\partial\chi(\bar{x}\mid c)=\partial^{\infty}\chi(\bar{x}\mid
C)$ is the normal cone to $C$ at $\bar{x}$. $\widehat{N}_{C}(\bar{x})\subseteq
N_{C}(\bar{x})$ always holds, and is an equality if $C$ is convex.

If $C$ is convex, $N_{C}(\bar{x})$ is the closed convex cone of elements
$v\in\mathds{R}^{n}$ such that $v^{T}(x-\bar{x})\leq 0$ for all $x\in C$ [22,
Theorem $6.9$].

For a set-valued mapping $S:X\rightarrow U$ where, for any $x\in X$,
$S(x)\subset U$, the graph of $S$ is the set
$\mathop{\textup{gph}}S:=\\{(x,u)\mid u\in S(x)\\}$.

For $\bar{x}\in\mathds{R}^{n}$, the limit superior of $S$ at $\bar{x}$ is
$\limsup_{x\rightarrow\bar{x}}S(x):=\\{u\mid\exists\
\\{x_{k}\\}\rightarrow\bar{x},\exists\ \\{u_{k}\\}\rightarrow u\text{ with
}u_{k}\in S(x_{k})\\}$, and the limit inferior of $S$ at $\bar{x}$ is
$\liminf_{x\rightarrow\bar{x}}S(x):=\\{u\mid\forall\
\\{x_{k}\\}\rightarrow\bar{x},\exists\
N\in\mathcal{N}_{\infty},\\{u_{k}\\}\underset{N}{\rightarrow}u\text{ with
}u_{k}\in S(x_{k})\\}$. $S(\bar{x})\subseteq\limsup_{x\rightarrow\bar{x}}S(x)$
and $S(\bar{x})\supseteq\liminf_{x\rightarrow\bar{x}}S(x)$ always hold.

The set-valued mapping $S$ is outer semicontinuous (osc) at $\bar{x}$ if
$\limsup_{x\rightarrow\bar{x}}S(x)\subseteq S(\bar{x})$, or, equivalently,
$\limsup_{x\rightarrow\bar{x}}S(x)=S(\bar{x})$. It is inner semicontinuous
(isc) at $\bar{x}$ if $\liminf_{x\rightarrow\bar{x}}S(x)\supseteq S(\bar{x})$,
or equivalently $\liminf_{x\rightarrow\bar{x}}S(x)=S(\bar{x})$ when $S$ is
closed-valued. If both conditions hold, $S$ is continuous at $\bar{x}$, i.e.,
$S(x)\rightarrow S(\bar{x})$ as $x\rightarrow\bar{x}$.

[22, Proposition $8.7$] For $\phi:\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
and $\bar{x}$ where $\phi$ is finite, $\partial\phi$ is osc at $\bar{x}$ with
respect to $\phi(x)\rightarrow\phi(\bar{x})$ when $x\rightarrow\bar{x}$, i.e.
for any $\\{x_{k}\\}\rightarrow\bar{x}$ with
$\\{\phi(x_{k})\\}\rightarrow\phi(\bar{x})$, there exists
$v_{k}\in\partial\phi(x_{k})$ for all $k$ such that
$\\{v_{k}\\}\rightarrow\bar{v}\in\partial\phi(\bar{x})$.

The graphical outer limit of a sequence of set-valued mappings $S_{k}$ is
defined by $(\mathop{\textup{g-lim sup}}_{k}S_{k})(x):=\\{u\mid\exists
N\in\mathcal{N}_{\infty}^{\\#},\\{x_{k}\\}\underset{N}{\rightarrow}x,\\{u_{k}\\}\underset{N}{\rightarrow}u,u_{k}\in
S_{k}(x_{k})\\}$. The graphical inner limit of a sequence of set-valued
mappings $S_{k}$ is defined by $(\mathop{\textup{g-lim
inf}}_{k}S_{k})(x):=\\{u\mid\exists
N\in\mathcal{N}_{\infty},\\{x_{k}\\}\underset{N}{\rightarrow}x,\\{u_{k}\\}\underset{N}{\rightarrow}u,u_{k}\in
S_{k}(x_{k})\\}$. If both limits agree, the graphical limit
$S=\mathop{\textup{g-lim}}S_{k}$ exists, so that we can also write
$S_{k}\overset{\textup{g}}{\longrightarrow}S$, and we have
$S_{k}\overset{\textup{g}}{\longrightarrow}S\iff\mathop{\textup{gph}}S_{k}\rightarrow\mathop{\textup{gph}}S$.

The epigraph of $\phi$ is the set
$\mathop{\textup{epi}}\phi=\\{(x,\alpha)\in\mathds{R}^{n}\times\mathds{R}\mid\alpha\geq\phi(x)\\}$.

We denote $\mathop{\textup{cl}}(\phi)$ the (lower) closure of $\phi$, i.e.,
the largest function less than $\phi$ that is lower semi-continuous. Its
epigraph is the closure of $\mathop{\textup{epi}}\phi$.

If $\phi_{k}:\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
for all $k\geq 0$, the lower and upper pointwise limits of $\\{\phi_{k}\\}$
are the functions $\mathop{\textup{p-lim inf}}_{k}\phi_{k}$ and
$\mathop{\textup{p-lim sup}}_{k}\phi_{k}:\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
defined for all $x\in\mathds{R}^{n}$ by

$\displaystyle(\mathop{\textup{p-lim inf}}_{k}\phi_{k})(x)$
$\displaystyle:=\liminf_{k}\phi_{k}(x),$ $\displaystyle(\mathop{\textup{p-lim
sup}}_{k}\phi_{k})(x)$ $\displaystyle:=\limsup_{k}\phi_{k}(x).$

When $\mathop{\textup{p-lim inf}}_{k}\phi_{k}$ and $\mathop{\textup{p-lim
sup}}_{k}\phi_{k}$ coincide, their common value is the pointwise limit of
$\\{\phi_{k}\\}$ denoted $\mathop{\textup{p-lim}}_{k}\phi_{k}$. If
$\mathop{\textup{p-lim}}_{k}\phi_{k}=\phi$, we write
$\\{\phi_{k}\\}\overset{\textup{p}}{\longrightarrow}\phi$.

For a sequence $\\{E_{k}\\}\subseteq\mathds{R}^{n}$, we define

$\displaystyle\limsup_{k}E_{k}$
$\displaystyle:=\\{z\in\mathds{R}^{n}\mid\exists
N\in\mathcal{N}_{\infty}^{\\#},\exists\\{z_{j}\\}_{N}\to z,\ z_{j}\in
E_{j}\text{ for all }j\in N\\}$ $\displaystyle\liminf_{k}E_{k}$
$\displaystyle:=\\{z\in\mathds{R}^{n}\mid\exists
N\in\mathcal{N}_{\infty},\exists\\{z_{j}\\}_{N}\to z,\ z_{j}\in E_{j}\text{
for all }j\in N\\}.$

Consider in particular $E_{k}:=\mathop{\textup{epi}}\phi_{k}$. It is not
difficult to see that $\limsup_{k}E_{k}$ and $\liminf_{k}E_{k}$ are also
epigraphs in the sense that if $(x,t)$ is in either set, then so is $(x,s)$
for any $s\geq t$. The lower epi-limit of $\\{\phi_{k}\\}$ is the function
$\mathop{\textup{e-lim inf}}_{k}\phi_{k}$ whose epigraph is
$\limsup_{k}E_{k}$, and the upper epi-limit of $\\{\phi_{k}\\}$ is the
function $\mathop{\textup{e-lim sup}}_{k}\phi_{k}$ whose epigraph is
$\liminf_{k}E_{k}$. It is always true that $\mathop{\textup{e-lim
inf}}_{k}\phi_{k}\leq\mathop{\textup{e-lim sup}}_{k}\phi_{k}$. When the two
coincide, their common value is called the epi-limit of $\\{\phi_{k}\\}$
denoted $\mathop{\textup{e-lim}}_{k}\phi_{k}$. If
$\mathop{\textup{e-lim}}_{k}\phi_{k}=\phi$, we also write
$\\{\phi_{k}\\}\overset{\textup{e}}{\longrightarrow}\phi$.

The following result summarizes important properties of epi-limits used in the
sequel.

[22, Proposition $7.4$] Let $\phi_{k}:\mathds{R}^{n}\to\mathds{R}$ for $k\geq
0$.

  1. 1.

$\mathop{\textup{e-lim inf}}_{k}\phi_{k}$ and $\mathop{\textup{e-lim
sup}}_{k}\phi_{k}$ are lsc, and so is $\mathop{\textup{e-lim}}_{k}\phi_{k}$
when it exists;

  2. 2.

if $\phi_{k}\geq\phi_{k+1}$ for all $k$, $\mathop{\textup{e-lim}}_{k}\phi_{k}$
exists and equals $\mathop{\textup{cl}}(\inf_{k}\phi_{k})$;

  3. 3.

if $\phi_{k}\leq\phi_{k+1}$ for al $k$, $\mathop{\textup{e-lim}}_{k}\phi_{k}$
exists and equals $\mathop{\textup{cl}}(\sup_{k}\phi_{k})$.

In addition, if $\phi$, $\underline{\phi}_{k}$, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{k}:\mathds{R}^{n}\to\mathds{R}$ with
$\underline{\phi}_{k}\leq\phi_{k}\leq\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{k}$ for all $k$, and if
$\\{\underline{\phi}_{k}\\}\overset{\textup{e}}{\longrightarrow}\phi$ and
$\\{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{k}\\}\overset{\textup{e}}{\longrightarrow}\phi$, then
$\\{\phi_{k}\\}\overset{\textup{e}}{\longrightarrow}\phi$.

The model that we will use in our algorithm uses an approximation
$\psi(\cdot;x)$ of $h$ at $x$ so that $\psi(s,x)\approx h(x+s)$. For
$\psi:\mathds{R}^{n}\times\mathds{R}^{m}\rightarrow\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$,
the function-valued mapping $x\mapsto\psi(\cdot;x)$ is _epi-continous_ at
$\bar{x}$ if
$\psi(\cdot,x)\overset{\textup{e}}{\longrightarrow}\psi(\cdot,\bar{x})$ as
$x\rightarrow\bar{x}$.

$\phi$ is level-bounded if, for every $\alpha\in\mathds{R}$, the lower level
set
$\mathop{\textup{lev}}_{\leq\alpha}\phi:=\\{x\in\mathds{R}^{n}\mid\phi(x)\leq\alpha\\}$
is bounded (possibly empty). The sequence of functions $\\{\phi_{k}\\}$ is
eventually level-bounded if, for each $\alpha\in\mathds{R}$, the sequence of
sets $\\{\mathop{\textup{lev}}_{\leq\alpha}\phi_{k}\\}$ is eventually bounded,
i.e., there is an index set $N\in\mathcal{N}_{\infty}$ such that
$\\{\mathop{\textup{lev}}_{\leq\alpha}\phi_{k}\\}_{k\in N}$ is bounded.

The following theorem establishes properties about the minimization of
sequences of epi-convergent functions.

[22, Theorem $7.33$] Suppose the sequence $\\{\phi_{k}\\}$ is eventually
level-bounded, and $\phi_{k}\overset{\textup{e}}{\longrightarrow}\phi$ with
$\phi_{k}$ and $\phi$ lsc and proper. Then,

(4) $\inf\phi_{k}\rightarrow\inf\phi$

with $-\infty<\inf\phi<+\infty$, while there exists $N\in\mathcal{N}_{\infty}$
such that $\mathop{\textup{argmin}}\phi_{k}$ is a bounded sequence of nonempty
sets with

(5)
$\limsup_{k}(\mathop{\textup{argmin}}\phi_{k})\subset\mathop{\textup{argmin}}\phi.$

Indeed, for any $\\{\epsilon_{k}\\}\searrow 0$ and
$x_{k}\in\epsilon_{k}\text{-}\mathop{\textup{argmin}}\phi_{k}$, $\\{x_{k}\\}$
is bounded and all its cluster points belong to
$\mathop{\textup{argmin}}\phi$. If $\mathop{\textup{argmin}}\phi$ consists of
a unique point $\bar{x}$, one must actually have
$\\{x_{k}\\}\rightarrow\bar{x}$.

The proximal operator associated with the proper lsc function $h$ and
parameter $\nu>0$ is

(6) $\mathop{\textup{prox}}_{\nu
h}(x):=\arg\min_{w}\tfrac{1}{2}\nu^{-1}\|w-x\|^{2}+h(w).$

##  3 Stationarity

First-order stationarity conditions for (1) may be stated as [22, Theorem
$10.1$]

(7) $0\in\nabla f(x)+\partial(h+\chi(\cdot\mid\mathds{R}^{n}_{+}))(x).$

We say that the constraint qualification (CQ) holds at $x$ if
$\partial^{\infty}(f+h)(x)$ contains no $v\neq 0$ such that $-v\in
N_{\mathds{R}^{n}_{+}}(x)$.

Under the constraint qualification, (7) can also be written [22, Theorem
$8.15$]

(8) $0\in\nabla f(x)+\partial h(x)+N_{\mathds{R}^{n}_{+}}(x).$

Our assumptions that $f$ is continuously differentiable and that $h$ is proper
lsc allows us to write [22, Exercise $8.8$ and Theorem $8.9$]

$\partial^{\infty}(f+h)(x)=\partial^{\infty}h(x)=\\{v\in\mathds{R}^{n}\mid(v,0)\in
N_{\mathop{\textup{epi}}h}(x,h(x))\\}.$

Due to the simple form of $N_{\mathds{R}^{n}_{+}}$, the constraint
qualification states that the only $v\in\partial^{\infty}h(x)$ such that
$v\geq 0$ and satisfying $v_{i}=0$ if $x_{i}>0$ is $v=0$.

A simple example where the constraint qualification is not satisfied is found
by setting $n=1$, $f=0$, and $h(x)=|x|_{0}$, i.e., $h(x)=1$ if $x\neq 0$ and
$h(0)=0$, in (1). The unique solution is $\bar{x}=0$. Then,
$\widehat{N}_{\mathop{\textup{epi}}h}(0,0)=N_{\mathop{\textup{epi}}h}(0,0)=\\{(v,t)\mid
t\leq 0\\}$ so that $\partial^{\infty}h(0)=\mathds{R}$. The qualification
condition requires that the only $v\in\mathds{R}$ such that $v\geq 0$ be
$v=0$, which is clearly not the case. Of course, the bound constraint in the
example above is redundant and the constrained and unconstrained solutions
coincide.

Using [22, Example $6.10$],
$N_{\mathds{R}_{+}^{n}}(x)=N_{\mathds{R}_{+}}(x_{1})\times\dots\times
N_{\mathds{R}_{+}}(x_{n})$, where $N_{\mathds{R}_{+}}(0)=(-\infty,0]$ and for
all $x_{i}>0$, $N_{\mathds{R}_{+}}(x_{i})=\\{0\\}$. Thus, (8) can also be
formulated as

(9) $0\in\nabla f(x)+\partial h(x)-z,\quad Xz=0,\quad x\geq 0,\quad z\geq 0,$

where $X=\mathop{\textup{diag}}(x)$ and $Z=\mathop{\textup{diag}}(z)$.

For fixed $x\in\mathds{R}^{n}_{+}$ and $z\in\mathds{R}_{+}^{n}$, we define
approximations

(10a) $\displaystyle\varphi^{\mathcal{L}}(s;x,z)$ $\displaystyle:=f(x)+(\nabla
f(x)-z)^{T}s,$ (10b) $\displaystyle\psi(s;x)$ $\displaystyle\phantom{:}\approx
h(x+s)\text{ with }\psi(0;x)=h(x)\text{ and }\partial\psi(0;x)=\partial h(x),$
(10c) $\displaystyle\hat{\psi}(s;x)$
$\displaystyle:=\psi(s;x)+\chi(x+s\mid\mathds{R}^{n}_{+}),$

and the model of $f+h$ about $x$

(11)
$m^{\mathcal{L}}(s;x,z,\nu):=\varphi^{\mathcal{L}}(s;x,z)+\tfrac{1}{2}\nu^{-1}\|s\|^{2}+\psi(s;x),$

where $\nu>0$. We point out that $\nabla\varphi^{\mathcal{L}}(s;x,z)=\nabla
f(x)-z$, which is the expression of the Lagrangian in the smooth case, thus,
we use the superscript$~{}^{\mathcal{L}}$ to denote objects sharing
similarities with the smooth Lagrangian.

For $\Delta\geq 0$, we further define

(12a) $\displaystyle p^{\mathcal{L}}(\Delta;x,z,\nu)$
$\displaystyle:=\min_{s}\
\varphi^{\mathcal{L}}(s;x,z)+\tfrac{1}{2}\nu^{-1}\|s\|^{2}+\hat{\psi}(s;x)+\chi(s\mid\Delta\mathds{B}),$
(12b) $\displaystyle P^{\mathcal{L}}(\Delta;x,z,\nu)$
$\displaystyle:=\mathop{\textup{argmin}}_{s}\
\varphi^{\mathcal{L}}(s;x,z)+\tfrac{1}{2}\nu^{-1}\|s\|^{2}+\hat{\psi}(s;x)+\chi(s\mid\Delta\mathds{B}).$

Our associated optimality measure is

(13)
$\xi^{\mathcal{L}}(\Delta;x,z,\nu):=f(x)+h(x)-\varphi^{\mathcal{L}}(s^{\mathcal{L}};x,z)-\psi(s^{\mathcal{L}};x),$

where $s^{\mathcal{L}}\in P^{\mathcal{L}}(\Delta,x,z,\nu)$.

Let the CQ hold at $x\in\mathds{R}^{n}_{+}$, $z\geq 0$ such that $Xz=0$ and
$\Delta>0$. Then, $\xi^{\mathcal{L}}(\Delta;x,z,\nu)=0$ $\Longleftrightarrow$
$0\in P^{\mathcal{L}}(\Delta;x,z,\nu)$ $\Longrightarrow$ $x$ is first-order
stationary for (1).

###### Proof.

The first equivalence follows directly from (12b)–(13). The first-order
necessary conditions for (12a) then imply

(14) $\displaystyle 0$
$\displaystyle\in\nabla\varphi^{\mathcal{L}}(0;x,z)+\partial(\hat{\psi}(\cdot;x)+\chi(\cdot\mid\Delta\mathds{B}))(0)$
$\displaystyle=\nabla
f(x)-z+\partial(\psi(\cdot;x)+\chi(\cdot\mid(-x+\mathds{R}_{+}^{n}))+\chi(\cdot\mid\Delta\mathds{B}))(0)$
$\displaystyle=\nabla
f(x)-z+\partial(\psi(\cdot;x)+\chi(\cdot\mid(-x+\mathds{R}_{+}^{n})\cap\Delta\mathds{B}))(0).$

As $(-x+\mathds{R}_{+}^{n})$ and $\Delta\mathds{B}$ are convex, so is
$(-x+\mathds{R}_{+}^{n})\cap\Delta\mathds{B}$. From this observation, we
deduce that,

$\displaystyle N_{(-x+\mathds{R}^{n}_{+})\cap\Delta\mathds{B}}(0)$
$\displaystyle=\partial\chi(0\mid(-x+\mathds{R}^{n}_{+})\cap\Delta\mathds{B})$
$\displaystyle=\partial\chi(0\mid(-x+\mathds{R}^{n}_{+}))+\partial\chi(0\mid\Delta\mathds{B})$
$\displaystyle=\partial\chi(x\mid\mathds{R}^{n}_{+})$
$\displaystyle=N_{\mathds{R}_{+}^{n}}(x).$

The CQ combined with the above equations indicate that there is no
$v\in\partial^{\infty}h(x)=\partial^{\infty}\psi(0;x)$, $v\neq 0$ such that
$-v\in
N_{\mathds{R}_{+}^{n}}(x)=N_{(-x+\mathds{R}_{+}^{n})\cap\Delta\mathds{B}}(0)$,
thus, [22, Corollary $10.9$] leads to

(15)
$\displaystyle\partial(\psi(\cdot;x)+\chi(\cdot\mid(-x+\mathds{R}_{+}^{n})\cap\Delta\mathds{B}))(0)$
$\displaystyle\subset\partial\psi(0;x)+N_{(-x+\mathds{R}_{+}^{n})\cap\Delta\mathds{B}}(0)$
$\displaystyle=\partial h(x)+N_{\mathds{R}^{n}_{+}}(x).$

By injecting (15) into (14), we obtain

$0\in\nabla f(x)-z+\partial h(x)+N_{\mathds{R}^{n}_{+}}(x).$

From the observation above (9) and the fact that $Xz=0$, we deduce that for
any $v\in N_{\mathds{R}^{n}_{+}}(x)$, $v-z\in N_{\mathds{R}_{+}^{n}}(x)$.
Thus,

$0\in\nabla f(x)+\partial h(x)+N_{\mathds{R}_{+}^{n}}(x).$

If $h$ is convex, the CQ is not required in Section 3 [22, Exercise $10.8$].

##  4 Projected-directions methods

Let us briefly recall the proximal gradient method [19] used to solve

(16) $\mathop{\textup{minimize}}_{s\in\mathds{R}^{n}}\ f(s)+\tilde{h}(s),$

where $f:\mathds{R}^{n}\to\mathds{R}$ has Lipschitz-continuous gradient and
$\tilde{h}:\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
is proper and lower semi-continuous. The method generates iterates $s_{k}$
such that

(17) $s_{k+1}\in\mathop{\textup{prox}}_{\nu\tilde{h}}(s_{k}-\nu\nabla
f(s_{k}))=\mathop{\textup{argmin}}_{s}f(s_{k})+\nabla
f(s_{k})^{T}(s-s_{k})+\tfrac{1}{2}\nu^{-1}\|s-s_{k}\|^{2}+\tilde{h}(s),$

where $\nu>0$, which leads to the first-order stationarity conditions

(18) $0\in s_{k+1}-s_{k}+\nu\nabla f(s_{k})+\nu\partial\tilde{h}(s_{k+1}).$

A first approach to solving (1), that we can reformulate as

$\mathop{\textup{minimize}}_{x\in\mathds{R}^{n}}\
f(x)+h(x)+\chi(x\mid\mathds{R}^{n}_{+}),$

is to use projected-directions methods. A simple example of such methods
consists in performing the identification
$\tilde{h}=h+\chi(\cdot\mid\mathds{R}_{+}^{n})$ in (16), and using the
proximal gradient method as in (17).

Aravkin et al.’s TR and R2 are other examples of algorithms that can solve (1)
with a similar strategy. Replacing $h$ by $\tilde{h}$ and $\psi(\cdot;x)$ by
$\tilde{\psi}(\cdot;x)=\psi(\cdot;x)+\chi(\cdot\mid-x+\mathds{R}_{+}^{n})$ in
all models of TR and R2 is sufficient to generalize these methods to (1), if a
solution of

(19)
$s_{k}\in\mathop{\textup{prox}}_{\nu_{k}\tilde{\psi}(\cdot;x_{k})+\chi(\cdot\mid\Delta_{k}\mathds{B})}(-\nu_{k}\nabla
f(x_{k})).$

for TR, or

(20)
$s_{k}\in\mathop{\textup{prox}}_{\nu_{k}\tilde{\psi}(\cdot;x_{k})}(-\nu_{k}\nabla
f(x_{k}))$

for R2, is available.

Leconte and Orban [17] implement a variant of TR named TRDH that handles bound
constraints for separable regularizers $h$ (assuming that $\psi$ is also
separable). TRDH solves at each iteration $k$ and for all
$i\in\\{1,\ldots,n\\}$ the problem

(21) $(s_{k})_{i}\in\arg\min_{s_{i}}\nabla
f(x_{k}))_{i}s_{i}+\tfrac{1}{2}(d_{k})_{i}s_{i}^{2}+(\psi(s;x_{k}))_{i}+\chi(s_{i}\mid\Delta_{k}\mathds{B}\cap(-(x_{k})_{i}+\mathds{R}_{+}),$

with $(d_{k})_{i}\in\mathds{R}$. The special choice $(d_{k})_{i}=\nu_{k}^{-1}$
shows that solving (21) for all $i$ is equivalent to solving (19).

However, for nonseperable regularizers $h$, projected-directions methods rely
on computing search directions such as (19), (20) or (21), which may be
complicated (impossible for the latter), and therefore seems to be a
limitation of this approach. The following section describes the
implementation of a method that is different from projected directions
methods, and is based upon interior-point techniques.

##  5 Barrier methods

Consider a sequence $\\{\mu_{k}\\}\searrow 0$.

Let $\phi_{k}$ be defined as in (3). Then,
$\mathop{\textup{e-lim}}\phi_{k}=\chi(\cdot\mid\mathds{R}^{n}_{+})$.

###### Proof.

It is sufficient to show that
$\mathop{\textup{e-lim}}\phi_{k,i}=\chi(\cdot\mid\mathds{R}_{+})$ for
$i=1,\ldots,n$. Our goal is to bound each $\phi_{k}$ by two functions having
$\chi(\cdot\mid\mathds{R}_{+}^{n})$ as epi-limit. We define

$\phi_{k,i}^{>}(x):=\begin{cases}+\infty&\text{if }x\leq 0\\\
\phi_{k,i}(x)&\text{if }0<x<1\\\ 0&\text{if }x\geq
1,\end{cases}\qquad\phi_{k,i}^{<}(x):=\begin{cases}+\infty&\text{if }x\leq
0\\\ 0&\text{if }0<x<1\\\ \phi_{k,i}(x)&\text{if }x\geq 1.\end{cases}$

By construction, $\phi_{k,i}(x)=\phi_{k,i}^{>}(x)+\phi_{k,i}^{<}(x)$,
$\\{\phi_{k,i}^{>}(x)\\}\searrow 0$ and $\\{\phi_{k,i}^{<}(x)\\}\uparrow 0$ as
$k\to\infty$ for all $x>0$. In particular,
$\\{\phi_{k,i}^{>}\\}\overset{\textup{p}}{\longrightarrow}\chi(\cdot\mid\mathds{R}_{+})$
and
$\\{\phi_{k,i}^{<}\\}\overset{\textup{p}}{\longrightarrow}\chi(\cdot\mid\mathds{R}_{+})$
as $k\to\infty$.

By [22, Proposition $7.4$c], because $\\{\phi_{k,i}^{>}\\}$ is nonincreasing
with $k$, its epi-limit is well defined and
$\\{\phi_{k,i}^{>}\\}\overset{\textup{e}}{\longrightarrow}\mathop{\textup{cl}}\inf_{k}\phi_{k,i}^{>}=\chi(\cdot\mid\mathds{R}_{+})$.

Similarly, by [22, Proposition $7.4$d], because $\\{\phi_{k,i}^{<}\\}$ is
nondecreasing with $k$, its epi-limit is well defined and
$\phi_{k,i}^{<}\overset{\textup{e}}{\longrightarrow}\sup_{k}\mathop{\textup{cl}}\phi_{k,i}^{<}=\chi(\cdot\mid\mathds{R}_{+})$.

Because $\phi_{k,i}^{<}\leq\phi_{k,i}\leq\phi_{k,i}^{>}$, [22, Proposition
$7.4$g] implies that
$\\{\phi_{k,i}\\}\overset{\textup{e}}{\longrightarrow}\chi(\cdot\mid\mathds{R}_{+})$,
and consequently, we obtain
$\\{\phi_{k}\\}\overset{\textup{e}}{\longrightarrow}\chi(\cdot\mid\mathds{R}^{n}_{+})$.

$\mathop{\textup{e-lim}}f+h+\phi_{k}=f+h+\chi(\cdot\mid\mathds{R}^{n}_{+})$.

###### Proof.

Section 5 and [22, Theorem $7.46$a] imply
$\\{h+\phi_{k}\\}\overset{\textup{e}}{\longrightarrow}h+\chi(\cdot\mid\mathds{R}^{n}_{+})$.
Finally, because $f$ is continuous, [22, Exercise $7.8$a] yields
$\\{f+h+\phi_{k}\\}\overset{\textup{e}}{\longrightarrow}f+h+\chi(\cdot\mid\mathds{R}^{n}_{+})$.

The following corollary legitimizes the barrier approach for (1).

Let $\inf\\{f(x)+h(x)\mid x\geq 0\\}$ be finite. For all $\epsilon\geq 0$,

$\limsup_{k}(\epsilon\text{-}\mathop{\textup{argmin}}f+h+\phi_{k})\subset\epsilon\text{-}\mathop{\textup{argmin}}f+h+\chi(\cdot\mid\mathds{R}^{n}_{+}).$

In particular, if $\\{\epsilon_{k}\\}\searrow 0$,

$\limsup_{k}(\epsilon_{k}\text{-}\mathop{\textup{argmin}}f+h+\phi_{k})\subset\mathop{\textup{argmin}}f+h+\chi(\cdot\mid\mathds{R}^{n}_{+}).$

###### Proof.

Follows directly from [22, Theorem $7.31$b].

If $\inf\\{f(x)+h(x)\mid x\geq 0\\}$ is finite, the definition of the limit
superior of a sequence of sets and the second part of Section 5 indicate that
for any $\\{\epsilon_{k}\\}\searrow 0$, there exists
$N\in\mathcal{N}_{\infty}^{\\#}$ and
$\bar{x}_{k}\in\epsilon_{k}\text{-}\mathop{\textup{argmin}}f+h+\phi_{k}$ for
all $k\in N$ such that $\\{\bar{x}_{k}\\}_{N}$ converges to a solution of (1).

###  5.1 Barrier subproblem

The $k$-th subproblem is

(22) $\mathop{\textup{minimize}}_{x}\ f(x)+h(x)+\phi_{k}(x).$

$x_{k}^{*}\in\mathds{R}^{n}_{++}$ is first-order stationary for (22) if

(23) $0\in\nabla f(x_{k}^{*})-\mu_{k}(X_{k}^{*})^{-1}e+\partial h(x_{k}^{*}).$

We call the process of solving (22) the $k$-th sequence of _inner_ iterations,
and we denote its iterates $x_{k,j}$ for $j\geq 0$. The definition of (22)
along with certain parameter updates will be called an _outer_ iteration.

For $x\in\mathds{R}^{n}_{++}$ and $\delta\in(0,1)$, let

(24)
$\mathds{R}^{n}_{\delta}(x):=\\{s\in\mathds{R}^{n}\mid\min_{i}(x+s)_{i}\geq\delta\min_{i}x_{i}\\}\subset(-x+\mathds{R}^{n}_{++}).$

Note that $\mathds{R}^{n}_{\delta}(x)$ is closed, and also convex, as shown in
the following lemma.

Let $\delta\in(0,1)$ and $x\in\mathds{R}_{++}^{n}$. Then
$\mathds{R}_{\delta}^{n}(x)$ defined in (24) is convex.

###### Proof.

Let $s_{1}$ and $s_{2}\in\mathds{R}_{\delta}^{n}(x)$, and $t\in[0,1]$. By
definition, $\min_{i}(s_{1}+x)_{i}\geq\delta\min_{i}x_{i}$ and
$\min_{i}(s_{2}+x)_{i}\geq\delta\min_{i}x_{i}$. Now,

$\displaystyle\min_{i}(ts_{1}+(1-t)s_{2}+x)_{i}$
$\displaystyle\geq\min_{i}(t(s_{1}+x))_{i}+\min_{i}((1-t)(s_{2}+x))_{i}$
$\displaystyle=t\min_{i}(s_{1}+x)_{i}+(1-t)\min_{i}(s_{2}+x)_{i}$
$\displaystyle\geq t\delta\min_{i}x_{i}+(1-t)\delta\min_{i}x_{i}$
$\displaystyle=\delta\min_{i}x_{i}.$

Thus, $ts_{1}+(1-t)s_{2}\in\mathds{R}_{\delta}^{n}(x)$ and
$\mathds{R}_{\delta}^{n}(x)$ is convex.

Under the assumptions of Section 5.1 and for $\Delta>0$,
$\Delta\mathds{B}\cap\mathds{R}_{\delta}^{n}(x)$ is convex.

At outer iteration $k$, we choose $\delta_{k}\in(0,1)$, and solve (22)
inexactly by approximately solving a sequence of trust-region subproblems of
the form

(25a) $\displaystyle\mathop{\textup{minimize}}_{s}$ $\displaystyle\
m(s;x_{k,j})+\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}^{n}_{\delta_{k}}(x_{k,j})),$
(25b) $\displaystyle m(s;x_{k,j})$
$\displaystyle:=\varphi(s;x_{k,j})+\psi(s;x_{k,j}),$

where $\varphi(s;x_{k,j})\approx(f+\phi_{k})(x_{k,j}+s)$ and
$\psi(s;x_{k,j})\approx h(x_{k,j}+s)$ model the smooth and nonsmooth parts of
(22), respectively, and $\Delta_{k,j}>0$ is a trust-region radius. Models are
required to satisfy the following assumption.

######  Model Assumption 5.1.

For any $k\in\mathds{N}$ and $j\in\mathds{N}$, $\varphi(\cdot;x_{k,j})$ is
continuously differentiable on $\mathds{R}^{n}_{++}$ with
$\varphi(0;x_{k,j})=f(x_{k,j})+\phi_{k}(x_{k,j})$ and
$\nabla_{s}\varphi(0;x_{k,j})=\nabla f(x_{k,j})+\nabla\phi_{k}(x_{k,j})$. In
addition, $\nabla_{s}\varphi(\cdot;x_{k,j})$ is Lipschitz continuous with
constant $L_{k,j}\geq 0$. We require that $\psi(\cdot;x_{k,j})$ be proper,
lsc, and satisfy $\psi(0;x_{k,j})=h(x_{k,j})$ and
$\partial\psi(0;x_{k,j})=\partial h(x_{k,j})$.

Let Model Assumption 5.1 be satisfied. Then $s=0$ is first-order stationary
for (25) if and only if $x_{k,j}$ is first-order stationary for (22).

###### Proof.

If $x_{k,j}$ is first-order stationary, then $0\in\nabla
f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e+\partial
h(x_{k,j})=\nabla_{s}\varphi(0;x_{k,j})+\partial\psi(0;x_{k,j})$. Note that
$x_{k,j}>0$ so
$\partial\chi(0\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}^{n}_{\delta_{k}}(x_{k,j}))=\\{0\\}$,
and $\Delta_{k,j}>0$, so there is an open set
$O\subset\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})$ such
that for all $s\in O$,
$\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}^{n}_{\delta_{k}}(x_{k,j}))=0$.
Thus, using the definition of the subdifferential,
$\partial(\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}^{n}_{\delta_{k}}(x_{k,j})))(0)=\partial\psi(0;x_{k,j})$.
We conclude that
$0\in\nabla_{s}\varphi(0;x_{k,j})+\partial(\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}^{n}_{\delta_{k}}(x_{k,j})))(0)$,
which is the definition of $s=0$ being first-order stationary for (25). The
reciprocal can also be established from these observations, because if
$0\in\nabla_{s}\varphi(0;x_{k,j})+\partial(\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}^{n}_{\delta_{k}}(x_{k,j})))(0)$,
we have shown that
$0\in\nabla_{s}\varphi(0;x_{k,j})+\partial\psi(0;x_{k,j})=\nabla
f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e+\partial h(x_{k,j})$.

As a special case of Section 5.1, if $s=0$ solves (25), then $x_{k,j}$ is
first-order stationary for (22).

Let

(26) $\varphi_{f}(s;x_{k,j},B_{k,j}):=f(x_{k,j})+\nabla
f(x_{k,j})^{T}s+\tfrac{1}{2}s^{T}B_{k,j}s,$

where $B_{k,j}=B_{k,j}^{T}$, be a second order Taylor approximation of $f$
about $x_{k,j}$. We are particularly interested in the quadratic model

(27) $\displaystyle\varphi(s;x_{k,j},B_{k,j})$
$\displaystyle:=\varphi_{f}(s;x_{k,j},B_{k,j})+\phi_{k}(x_{k,j})-\mu_{k}e^{T}X_{k,j}^{-1}s+\tfrac{1}{2}s^{T}X_{k,j}^{-1}Z_{k,j}s$
$\displaystyle=(f+\phi_{k})(x_{k,j})+(\nabla
f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e)^{T}s+\tfrac{1}{2}s^{T}(B_{k,j}+X_{k,j}^{-1}Z_{k,j})s,$

where $z_{k,j}$ is an approximation to the vector of multipliers for the bound
constraints of (1).

Let $s_{k,j}$ be an approximate solution of (25). If $s_{k,j}$ is accepted as
a step for our algorithm used to solve (22) (the acceptance condition is
detailed in Algorithm 2), we perform the update $x_{k,j+1}=x_{k,j}+s_{k,j}$.

By analogy with the smooth case, we use $z_{k,j}:=\mu_{k}X_{k,j}^{-1}e$ when
$x_{k,j}$ is first-order stationary for (22). Multiplying through by
$X_{k,j}$, we obtain $X_{k,j}z_{k,j}=\mu_{k}e$. Linearizing the continuous
equality $Xz-\mu e=0$ with respect to $x$ and $z$ and evaluating all
quantities at iteration $(k,j)$ yields $X_{k,j}\Delta
z_{k,j}+Z_{k,j}s_{k,j}=\mu_{k}e-X_{k,j}z_{k,j}$, which suggests that if
$x_{k,j+1}=x_{k,j}+s_{k,j}$, then $z_{k,j+1}=z_{k,j}+\Delta
z_{k,j}=\mu_{k}X_{k,j}^{-1}e-X_{k,j}^{-1}Z_{k,j}s_{k,j}$.

However, the latter $z_{k,j+1}$ may not be positive. We perform the update
described by Conn et al. [9], by defining

(28) $\hat{z}_{k,j+1}=\mu_{k}X_{k,j}^{-1}e-X_{k,j}^{-1}Z_{k,j}s_{k,j},$

and projecting $\hat{z}_{k,j+1}$ componentwise into the following interval to
get $z_{k,j+1}$

(29)
$\mathcal{I}=[\kappa_{\textup{zul}}\min(e,z_{k,j},\mu_{k}X_{k,j+1}^{-1}e),\
\max(\kappa_{\textup{zuu}}e,z_{k,j},\kappa_{\textup{zuu}}\mu_{k}^{-1}e,\kappa_{\textup{zuu}}\mu_{k}X_{k,j+1}^{-1}e)],$

with $0<\kappa_{\textup{zul}}<1<\kappa_{\textup{zuu}}$. Projecting
$\hat{z}_{k,j+1}$ into (29) always generates a positive $z_{k,j+1}$. The
choice $z_{k,j+1}=z_{k,j}$ is also available. The other bounds of (29) will be
useful in Section 5.2 and Section 5.3.

We define the following model, based upon a first-order Taylor approximation

(30a) $\displaystyle\varphi_{\textup{cp}}(s;x_{k,j})$
$\displaystyle:=(f+\phi_{k})(x_{k,j})+(\nabla
f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e)^{T}s,$ (30b) $\displaystyle
m_{\textup{cp}}(s;x_{k,j},\nu_{k,j})$
$\displaystyle:=\varphi_{\textup{cp}}(s;x_{k,j})+\tfrac{1}{2}\nu_{k,j}^{-1}\|s\|^{2}+\psi(s;x_{k,j}),$

where “cp” stands for “Cauchy point”. Let $s_{k,j,1}$ be the solution of (25)
with model $m_{\textup{cp}}(s;x_{k,j},\nu_{k,j})$. As stated in [2, Section
$3.2$], $s_{k,j,1}$ is actually the first step of the proximal gradient method
(17) from $s_{k,j,0}=0$ applied to the minimization of
$\varphi_{\textup{cp}}+\psi$ with step length $\nu_{k,j}$:

(31)
$s_{k,j,1}\in\mathop{\textup{prox}}_{\nu_{k,j}\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))}(-\nu_{k,j}\nabla\varphi_{\textup{cp}}(0;x_{k,j},\nu_{k,j})).$

Let

(32)
$\xi_{\textup{cp}}(\Delta;x_{k,j},\nu_{k,j}):=(f+\phi_{k}+h)(x_{k,j})-(\varphi_{\textup{cp}}+\psi)(s_{k,j,1};x_{k,j}),$

where $\Delta>0$, and let
$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta;x_{k,j},\nu_{k,j})^{1/2}$ be our
measure of criticality. Aravkin et al. [2] and its corrigendum [1] indicate
that $\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta;x_{k,j},\nu_{k,j})^{1/2}$ is
similar to $\nu_{k,j}^{-1}\|s_{k,j,1}\|$, which is the norm of the generalized
gradient at $x_{k,j}$. We can apply [22, Theorems $1.17$ and $7.41$] to
conclude that $\xi_{\textup{cp}}(\Delta;x_{k,j},\nu_{k,j})$ is proper lsc in
$(x_{k,j},\nu_{k,j})\in\mathds{R}_{++}^{n}\times\mathds{R}_{++}$. In
particular,
$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta;x_{k,j},\nu_{k,j})^{1/2}=0$ for any
$\Delta>0$ and $\nu_{k,j}>0$ $\Longrightarrow$ $s=0$ solves (25), and
$x_{k,j}$ is first-order stationary for (22).

Algorithm 1 summarizes the outer iteration.

Algorithm 1 Nonsmooth interior-point method (outer iteration).

1:Choose $\epsilon>0$, sequences $\\{\mu_{k}\\}\searrow 0$,
$\\{\epsilon_{d,k}\\}\searrow 0$, $\\{\epsilon_{p,k}\\}\searrow 0$, and
$\\{\delta_{k}\\}\rightarrow\bar{\delta}\in[0,1)$ with $\delta_{k}\in(0,1)$
for all $k$.

2:Choose $x_{0,0}\in\mathds{R}^{n}_{++}$ where $h$ is finite.

3:for $k=0,1,\ldots$ do

4: Compute an approximate solution $x_{k}:=x_{k,j}$ to (22) and
$z_{k}:=z_{k,j}$ in the sense that  (33)
$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}\leq\epsilon_{d,k}$
and  (34) $\|X_{k,j}z_{k,j}-\mu_{k}e\|\leq\epsilon_{p,k}.$

5: Set $x_{k+1,0}:=x_{k}$.

6:end for

For each outer iteration $k$, the inner iterations generate a sequence
$\\{x_{k,j}\\}$ according to an adaptation of [2, Algorithm $3.1$] in which
the subproblems have the form (25) with the smooth part of the model defined
by (27). Each trust-region step is required to satisfy the following
assumption.

######  Step Assumption 5.1.

Let $k\in\mathds{N}$. There exists $\kappa_{\textup{m},k}>0$ and
$\kappa_{\textup{mdc},k}\in(0,\,1)$ such that for all $j$,
$s_{k,j}\in\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}(x_{k,j})$,

(35a)
$\displaystyle|(f+\phi_{k}+h)(x_{k,j}+s_{k,j})-m(s_{k,j};x_{k,j},B_{k,j})|$
$\displaystyle\leq\kappa_{\textup{m},k}\|s_{k,j}\|^{2},$ (35b) $\displaystyle
m(0;x_{k,j},B_{k,j})-m(s_{k,j};x_{k,j},B_{k,j})$
$\displaystyle\geq\kappa_{\textup{mdc},k}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j}),$

where $m$ is defined in (25b)–(27), and
$\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})$ is defined in (32).

In Step Assumption 5.1, the subscript “m” of $\kappa_{\textup{m},k}$ refers to
the model adequacy, and the subscript “mdc” of $\kappa_{\textup{mdc},k}$
refers to the model decrease.

Algorithm 2 summarizes the process.

Algorithm 2 Nonsmooth interior-point method (inner iteration).

1:Choose constants  $0<\eta_{1}\leq\eta_{2}<1,\hskip 433.62pt minus
433.62pt0<\gamma_{1}\leq\gamma_{2}<1<\gamma_{3}\leq\gamma_{4},\hskip 433.62pt
minus 433.62pt\Delta_{\max}>\Delta_{k,0}>0,\hskip 433.62pt minus
433.62pt\alpha>0,\,\beta\geq 1.$

2:Compute $f(x_{k,0})+\phi_{k}(x_{k,0})+h(x_{k,0})$.

3:for $j=0,1,\ldots$ do

4: Choose $0<\nu_{k,j}\leq 1/(L_{k,j}+\alpha^{-1}\Delta_{k,j}^{-1})$.

5: Define $m(s;x_{k,j},B_{k,j})$ as in (25b) satisfying Model Assumption 5.1.

6: Define $m_{\textup{cp}}(s;x_{k,j},\nu_{k,j})$ as in (30b).

7: Compute $s_{k,j,1}$ as a solution of (31).

8: Compute an approximate solution $s_{k,j}$ of (25) such that
$\|s_{k,j}\|\leq\min(\Delta_{k,j},\,\beta\|s_{k,j,1}\|)$.

9: Compute the ratio
$\rho_{k,j}:=\frac{(f+\phi_{k}+h)(x_{k,j})-(f+\phi_{k}+h)(x_{k,j}+s_{k,j})}{m(0;x_{k,j},B_{k,j})-m(s_{k,j};x_{k,j},B_{k,j})}.$

10: If $\rho_{k,j}\geq\eta_{1}$, set $x_{k,j+1}=x_{k,j}+s_{k,j}$ and update
$z_{k,j+1}$ according to (28) and (29). Otherwise, set $x_{k,j+1}=x_{k,j}$ and
$z_{k,j+1}=z_{k,j}$.

11: Update the trust-region radius according to
$\hat{\Delta}_{k,j+1}\in\left\\{\begin{array}[]{lll}{[\gamma_{3}\Delta_{k,j},\,\gamma_{4}\Delta_{k,j}]}&\text{
if }\rho_{k,j}\geq\eta_{2},&\text{(very successful iteration)}\\\
{[\gamma_{2}\Delta_{k,j},\,\Delta_{k,j}]}&\text{ if
}\eta_{1}\leq\rho_{k,j}<\eta_{2},&\text{(successful iteration)}\\\
{[\gamma_{1}\Delta_{k,j},\,\gamma_{2}\Delta_{k,j}]}&\text{ if
}\rho_{k,j}<\eta_{1}&\text{(unsuccessful iteration)}\end{array}\right.$ and
$\Delta_{k,j+1}=\min(\hat{\Delta}_{k,j+1},\,\Delta_{\max})$.

12:end for

###  5.2 Convergence of the inner iterations

Let $k$ and $j$ be fixed positive integers, and $x_{k,j}>0$. We may rewrite
our “Cauchy point” subproblem for the inner iterations as in [2]:

(36a) $\displaystyle p(\Delta;x_{k,j},\nu_{k,j},\delta_{k})$
$\displaystyle:=\mathop{\textup{minimize}}_{s}m_{\textup{cp}}(s;x_{k,j},\nu_{k,j})+\chi(s;\Delta\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))$
(36b) $\displaystyle P(\Delta;x_{k,j},\nu_{k,j},\delta_{k})$
$\displaystyle:=\mathop{\textup{argmin}}_{s}m_{\textup{cp}}(s;x_{k,j},\nu_{k,j})+\chi(s;\Delta\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})).$

First, we present some properties of the subproblem (36a) in the following
result.

[2, Proposition $3.1$] Let Model Assumption 5.1 be satisfied, $\nu>0$,
$\delta>0$ and $x\in\mathds{R}_{+}^{n}$. If we define
$p(0;x,\nu,\delta):=\varphi_{\textup{cp}}(0;x)+\psi(0;x)$ and
$P(0;x,\nu,\delta)=\\{0\\}$, the domain of $p(\cdot;x,\nu,\delta)$ and
$P(\cdot;x,\nu,\delta)$ is $\\{\Delta\mid\Delta\geq 0\\}$. In addition,

  1. 1.

$p(\cdot;x,\nu,\delta)$ is proper lsc and for each $\Delta\geq 0$,
$P(\Delta;x,\nu,\delta)$ is nonempty and compact;

  2. 2.

if $\\{\Delta_{k,j}\\}\rightarrow\bar{\Delta}_{k}\geq 0$ in such a way that
$\\{p(\Delta_{k,j};x,\nu,\delta)\\}\rightarrow
p(\bar{\Delta}_{k};x,\nu,\delta)$, and for each $j$, $s_{k,j}\in
P(\Delta_{k,j};x,\nu,\delta)$, then $\\{s_{k,j}\\}$ is bounded and all its
limit points are in $P(\bar{\Delta}_{k};x,\nu,\delta)$;

  3. 3.

if
$\varphi_{\textup{cp}}(\cdot;x)+\tfrac{1}{2}\nu^{-1}\|s\|^{2}+\psi(\cdot;x)$
is strictly convex, $P(\Delta;x,\nu,\delta)$ is single valued;

  4. 4.

if $\bar{\Delta}_{k}>0$ and there exists $\bar{s}\in
P(\bar{\Delta}_{k};x,\nu,\delta)$ such that
$\bar{s}\in\mathrm{int}(\bar{\Delta}_{k}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x))$,
then $p(\cdot;x,\nu,\delta)$ is continuous at $\bar{\Delta}_{k}$ and
$\\{p(\Delta_{k,j};x,\nu,\delta)\\}\rightarrow
p(\bar{\Delta}_{k};x,\nu,\delta)$ holds in part 2.

###### Proof.

Model Assumption 5.1 and the compactness of
$\Delta\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x)$ ensure that the objective
of (36a) is always level-bounded in $s$ locally uniformly in $\Delta$, because
for any $\bar{\Delta}>0$, $\epsilon>0$, and
$\Delta\in(\bar{\Delta}-\epsilon,\bar{\Delta}+\epsilon)$ with $\Delta\geq 0$,
its level sets are contained in
$\Delta\mathds{B}\cap\mathds{R}_{\delta_{k}}(x_{k,j})\subseteq(\bar{\Delta}+\epsilon)\mathds{B}\cap\mathds{R}_{\delta_{k}}(x_{k,j})$.
From this observation, we can draw similar conclusions to the analysis of [2,
Proposition $3.1$].

The observation in the proof of Section 5.2 and Model Assumption 5.1 allows us
to derive directly some of the convergence properties of Aravkin et al. [2]
for Algorithm 2.

[2, Theorem $3.4$] Let Model Assumption 5.1 and Step Assumption 5.1 be
satisfied and let

(37)
$\Delta_{\textup{succ},k}:=\frac{\kappa_{\textup{mdc},k}(1-\eta_{2})}{2\kappa_{\textup{m},k}\alpha\beta^{2}}>0.$

If $x_{k,j}$ is not first-order stationary for (22) and
$\Delta_{k,j}\leq\Delta_{\textup{succ}}$, then iteration $j$ is very
successful and $\Delta_{k,j+1}\geq\Delta_{k,j}$.

###### Proof.

If $x_{k,j}$ is not first-order stationary,
$\mathds{R}_{\delta_{k}}^{n}(x_{k,j})\neq\\{0\\}$, thus $s_{k,j,1}\neq 0$ and
$s_{k,j}\neq 0$. The rest of the proof is identical to that of [2, Theorem
$3.4$].

Now, let
$\Delta_{\min,k}:=\min(\Delta_{k,0},\gamma_{1}\Delta_{\textup{succ},k})>0$.
Then, $\Delta_{k,j}\geq\Delta_{\min,k}$ for all $j\in\mathds{N}$. If we
consider $\varphi$ defined in (27), for $s_{1}$ and $s_{2}$ in
$\mathds{R}_{+}^{n}$,

(38)
$\|\nabla\varphi(s_{1})-\nabla\varphi(s_{2})\|=\|B_{k,j}(s_{1}-s_{2})+X_{k,j}^{-1}Z_{k,j}(s_{1}-s_{2})\|.$

As
$L_{k,j}=\|B_{k,j}+X_{k,j}^{-1}Z_{k,j}\|\leq\|B_{k,j}\|+\|X_{k,j}^{-1}\|\|Z_{k,j}\|$,
if $\\{B_{k,j}\\}_{j}$ remains bounded, $\\{x_{k,j}\\}_{j}$ must be bounded
away from zero to guarantee the existence of some $M_{k}>0$ such that
$L_{k,j}\leq M_{k}$ for all $j\in\mathds{N}$. To apply the complexity results,
and to establish that
$\liminf\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}=0$
if $f+h+\phi_{k}$ is bounded below on $\mathds{R}_{+}^{n}$, we need a stronger
assumption on the Lipschitz constant of our model.

######  Model Assumption 5.2 (2, Model Assumption $3.3$).

In Model Assumption 5.1, there exists $M_{k}>0$ such that $0\leq L_{k,j}\leq
M_{k}$ for all $j\in\mathds{N}$. In addition, we select $\nu_{k,j}$ in line 4
of Algorithm 2 in a way that there exists $\nu_{\min,k}>0$ such that
$\nu_{k,j}\geq\nu_{\min,k}$ for all $j\in\mathds{N}$.

As in Aravkin et al. [2], we can set
$\nu_{k,j}:=1/(L_{k,j}+\alpha^{-1}\Delta_{k,j}^{-1})$ in Algorithm 2 to ensure
$\nu_{k,j}\geq\nu_{\min,k}:=1/(M_{k}+\alpha^{-1}\Delta_{\min,k}^{-1})>0$ if
the first part of Model Assumption 5.2 holds with $M_{k}$. The observations
below (38) motivate us to prove that $\\{x_{k,j}\\}_{j}$ is bounded away from
zero in the next result.

Let $k\in\mathds{N}$, Model Assumption 5.1 be satisfied for $\varphi$ in (27),
and $(f+h)(x_{k,j})\geq(f+h)_{\textup{low},k}$. Then, there exists
$\kappa_{\textup{mdb},k}>0$ such that, for all $j$, we have

(39) $\min_{i}(x_{k,j})_{i}\geq\kappa_{\textup{mdb},k}.$

###### Proof.

We proceed similarly as in Conn et al. [9, Theorem $13.2.1$]. Let $k$ be a
positive integer and $\\{x_{k,j}\\}$ be a sequence generated by Algorithm 2.
As $\\{(f+\phi_{k}+h)(x_{k,j})\\}_{j}$ is decreasing and
$(f+h)(x_{k,j})\geq(f+h)_{\textup{low},k}$, we have
$\limsup_{j}\phi_{k}(x_{k,j})<\infty$, which implies that (39) holds.

The following proposition shows that we can use the convergence results of
Aravkin et al. [2] using the same assumptions they used for $\varphi_{f}$. It
will justify that Model Assumption 5.2 can be used for $\varphi$ defined in
(27).

Under the assumptions of Section 5.2, let $\varphi_{f}$ be defined as in (26)
so that $\nabla_{s}\varphi_{f}(\cdot;x_{k,j},B_{k,j})$ is Lipschitz continuous
with constant $\tilde{L}_{k,j}\geq 0$ and there exists $\tilde{M}_{k}>0$ such
that $0\leq\tilde{L}_{k,j}\leq\tilde{M}_{k}$ for all $j\in\mathds{N}$. Then
$\varphi$ satisfies Model Assumption 5.2.

###### Proof.

We can use (29) and (39) to say that $X_{k,j}^{-1}Z_{k,j}$ is bounded for all
$j$, and we deduce from (38) that $\varphi$ satisfies Model Assumption 5.2.

Now, we justify that Step Assumption 5.1 holds when
$h(x_{k,j}+s_{k,j})=\psi(s_{k,j};x_{k,j})$. As $\\{x_{k,j}\\}_{j}$ remains
bounded away from $\partial\mathds{R}_{+}^{n}$ with Section 5.2, so does
$\\{x_{k,j}+s_{k,j}\\}_{j}$ by definition of
$\mathds{R}_{\delta_{k}}(x_{k,j})$. Since $\nabla f$ is Lipschitz-continuous,

(40) $f(x_{k,j}+s_{k,j})-f(x_{k,j})-\nabla
f(x_{k,j})^{T}s_{k,j}\leq\tfrac{1}{2}L_{f}\|s_{k,j}\|^{2},$

and a second-order Taylor approximation of $\phi_{k}$ about $x_{k,j}+s_{k,j}$
gives

(41)
$\phi_{k}(x_{k,j}+s_{k,j})-\phi_{k}(x_{k,j})-\mu_{k}s_{k,j}^{T}X_{k,j}^{-1}e=\tfrac{1}{2}\mu_{k}s_{k,j}^{T}X_{k,j}^{-2}s_{k,j}+o(\|s_{k,j}\|^{2}).$

Under the assumptions of Section 5.2,
$L_{k,j}=\|B_{k,j}+X_{k,j}^{-1}Z_{k,j}\|\leq M_{k}$, and

$|(f+\phi_{k}+h)(x_{k,j}+s_{k,j})-m(s_{k,j};x_{k,j},B_{k,j})|\\\
=|f(x_{k,j}+s_{k,j})-f(x_{k,j})-\nabla
f(x_{k,j})^{T}s_{k,j}+\phi_{k}(x_{k,j}+s_{k,j})-\phi_{k}(x_{k,j})-\\\
\mu_{k}s_{k,j}^{T}X_{k,j}^{-1}e-\tfrac{1}{2}s_{k,j}^{T}(B_{k,j}+X_{k,j}^{-1}Z_{k,j})s_{k,j}|+o(\|s_{k,j}\|^{2}).$

The above equality combined with (40) and (41) implies that (35a) holds.

To emphasize the similarities between our inner iterations and the trust-
region algorithm of Aravkin et al. [2], and in light of Section 5.2, we use in
our next results that $\varphi$ satisfies Model Assumption 5.2, instead of
writting assumptions on $\varphi_{f}$. The following proposition gives us a
sufficient condition for (35b) to be satisfied.

[1, Proposition $1$] If Model Assumption 5.2 is satisfied with bounded Hessian
approximations $\\{B_{k,j}\\}_{j}$, then there exists
$\kappa_{\textup{mdc},k}\in(0,1)$ such that (35b) holds for all $j$.

###### Proof.

The proof is identical to that of [1, Proposition $1$] when replacing $B_{k}$
by $B_{k,j}+X_{k,j}^{-1}Z_{k,j}$, and using the subscripts$~{}_{k,j}$ where
$j$ is the iteration number of the algorithm instead of the
subscript$~{}_{k}$.

Section 5.2 allows us to write the following convergence results for Algorithm
2. As in Aravkin et al. [2], we define the smallest iteration number
$j_{k}(\epsilon)$ such that

(42)
$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}\leq\epsilon\quad(0<\epsilon<1),\quad
k=0,1,\ldots,$

and we express the set of successful iterations, the set of successful
iterations for which (42) has not yet been attained, and the set of
unsuccessful iterations for which (42) has not yet been attained as

(43a) $\displaystyle\mathcal{S}_{k}$
$\displaystyle:=\\{j\in\mathds{N}\mid\rho_{k,j}\geq\eta_{1}\\}$ (43b)
$\displaystyle\mathcal{S}_{k}(\epsilon)$
$\displaystyle:=\\{j\in\mathcal{S}_{k}\mid j<j_{k}(\epsilon)\\}$ (43c)
$\displaystyle\mathcal{U}_{k}(\epsilon)$ $\displaystyle:=\\{j\in\mathds{N}\mid
j\notin\mathcal{S}_{k}\text{ and }j<j_{k}(\epsilon)\\}.$

[2, Theorem $3.5$] Let Model Assumption 5.1 and Step Assumption 5.1 be
satisfied. If Algorithm 2 only generates finitely many successful iterations,
then $x_{k,j}=x_{k}^{*}$ for all sufficiently large $j$ and $x_{k}^{*}$ is
first-order critical for (22).

###### Proof.

The proof is inspired from [2, Theorem $3.5$], which itself follows that of
Conn et al. [9, Theorem $6.6.4$]. The assumptions indicate that there is
$j_{0}\in\mathds{N}$ such that all iterations $j\geq j_{0}$ are unsuccessful,
and $x_{k,j}=x_{k,j_{0}}=x_{k}^{*}$ because of the update rules of Algorithm
2. We assume by contradiction that $x_{k}^{*}$ is not first-order critical.
$x_{k}^{*}$ does not have any of its components equal to $+\infty$ because it
is attained after a finite number of iterations of Algorithm 2. As $h$ is
proper, $(f+h)(x_{k}^{*})>-\infty$. Thus, Section 5.2 implies that there
exists $\Delta_{k}^{*}$ such that
$\Delta_{k}^{*}\mathds{B}\subset\mathds{R}_{\delta_{k}}(x_{k}^{*})$. Since all
iterations $j\geq j_{0}$ are unsuccessful, there will be some $j_{1}\geq
j_{0}$ such that
$\Delta_{j_{1}}\leq\min(\Delta_{\textup{succ},k},\,\Delta_{k}^{*})$, which
implies that iteration $j_{1}$ is very successful with Section 5.2, and
contradicts the fact that $x_{k}^{*}$ is not first-order critical.

Finally, we have the following result for the inner iterates. [2, Theorem
$3.11$] Let $k\in\mathds{N}$, Model Assumption 5.1 and Model Assumption 5.2 be
verified for $\varphi$ defined in (27), and Step Assumption 5.1 be satisfied
for Algorithm 2. If there are infinitely many successful iterations, then,
either

(44)
$\underset{j\rightarrow\infty}{\lim}(f+\phi_{k}+h)(x_{k,j})=-\infty\quad\text{or}\quad\underset{j\rightarrow\infty}{\lim}\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}=0.$

###### Proof.

The proof is identical as that of [2, Theorem $3.11$].

Now, the definitions of $\varphi_{\textup{cp}}$ in (30a) and $s_{k,j,1}$ as a
minimizer of (25) with model $m_{\textup{cp}}$ indicate that

$(f+\phi_{k}+h)(x_{k,j})=\varphi_{\textup{cp}}(0;x_{k,j})+\psi(0;x_{k,j})=m_{\textup{cp}}(0;x_{k,j},\nu_{k,j})\geq\\\
m_{\textup{cp}}(s_{k,j,1};x_{k,j},\nu_{k,j})=\varphi_{\textup{cp}}(s_{k,j,1};x_{k,j})+\tfrac{1}{2}\nu_{k,j}^{-1}\|s_{k,j,1}\|^{2}+\psi(s_{k,j,1};x_{k,j}).$

By reinjecting this inequality into the definition of $\xi_{\textup{cp}}$ in
(32), we obtain

$\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})\geq\tfrac{1}{2}\nu_{k,j}^{-1}\|s_{k,j,1}\|^{2},$

so that

(45)
$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}\geq\tfrac{1}{\sqrt{2}}\nu_{k,j}^{-1}\|s_{k,j,1}\|.$

From this observation, we deduce that
$s_{k,j}\underset{j\rightarrow\infty}{\longrightarrow}0$ in the following
result. If Model Assumption 5.1 holds for $\varphi$ defined in (27) and

$\lim_{j\rightarrow\infty}\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}\rightarrow
0,$

then
$\lim_{j\rightarrow\infty}\|s_{k,j,1}\|=\lim_{j\rightarrow\infty}\|s_{k,j}\|=0$.

###### Proof.

We use $\beta\|s_{k,j,1}\|^{2}\geq\|s_{k,j}\|^{2}$ and (45) to conclude that

$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}\geq\tfrac{1}{\sqrt{2}}\nu_{k,j}^{-1}\|s_{k,j,1}\|\geq\tfrac{1}{\sqrt{2}}\nu_{k,j}^{-1}\beta^{-1/2}\|s_{k,j}\|\geq
0.$

With Model Assumption 5.1,
$\nu_{k,j}\underset{j\rightarrow\infty}{\longrightarrow}\bar{\nu}_{k}>0$, and
we have $\|s_{k,j}\|\underset{j\rightarrow\infty}{\longrightarrow}0$.

Now, we study the asymptotic satisfaction of the inner perturbed
complementarity. We show that (34) is eventually satisfied, similarly to [9,
Theorem $13.6.4$] in the smooth case.

Let $k\in\mathds{N}$, Model Assumption 5.1 and Model Assumption 5.2 be
verified for $\varphi$ defined in (27), Step Assumption 5.1 be satisfied for
Algorithm 2, and $(f+\phi_{k}+h)(x_{k,j})\geq(f+\phi_{k}+h)_{\textup{low},k}$
for all $j\in\mathds{N}$. Then,

$\lim_{j\rightarrow\infty}\|\mu_{k}X_{k,j}^{-1}e-z_{k,j}\|=0.$

###### Proof.

We proceed similarly as in the proof of [9, Theorem $13.6.4$]. With the
formula $\hat{z}_{k,j+1}=\mu_{k}X_{k,j}^{-1}e-X_{k,j}^{-1}Z_{k,j}s_{k,j}$,

$\displaystyle\|\hat{z}_{k,j+1}-\mu_{k}X_{k,j+1}^{-1}e\|$
$\displaystyle\leq\|X_{k,j}^{-1}Z_{k,j}s_{k,j}\|+\mu_{k}\|X_{k,j+1}^{-1}e-X_{k,j}^{-1}e\|$
$\displaystyle\leq\|X_{k,j}^{-1}Z_{k,j}\|\|s_{k,j}\|+\mu_{k}\sqrt{n}\|X_{k,j+1}^{-1}-X_{k,j}^{-1}\|.$

Using Model Assumption 5.2, we have
$\lim_{j\rightarrow\infty}\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}=0$.
This leads to $\lim_{j\rightarrow\infty}\|s_{k,j}\|=0$ with Equation 45, and
we have that $x_{k,j}$ is bounded away from $0$ using Section 5.2. Therefore,
either iteration $j$ is not successful and $X_{k,j+1}=X_{k,j}$, or

$\displaystyle\|X_{k,j+1}^{-1}-X_{k,j}^{-1}\|$
$\displaystyle=\|X_{k,j}^{-1}(X_{k,j}X_{k,j+1}^{-1}-I)\|$
$\displaystyle=\|X_{k,j}^{-1}(X_{k,j}(X_{k,j}+S_{k,j})^{-1}-I)\|$
$\displaystyle=\|X_{k,j}^{-1}((I+X_{k,j}^{-1}S_{k,j})^{-1}-I)\|\rightarrow 0.$

Since $X_{k,j+1}^{-1}z_{k,j}$ is also bounded for all $j$ using Section 5.2
and (29), we have

(46)
$\lim_{j\rightarrow\infty}\|\hat{z}_{k,j+1}-\mu_{k}X_{k,j+1}^{-1}e\|\rightarrow
0.$

For $j$ large enough, we have

(47)
$\kappa_{\textup{zul}}\mu_{k}X_{k,j+1}^{-1}e\leq\hat{z}_{k,j+1}\leq\kappa_{\textup{zuu}}\mu_{k}X_{k,j+1}^{-1}e,$

so that $\hat{z}_{k,j+1}=z_{k,j+1}$ if $j$ is large enough.

In Algorithm 1, each subproblem is solved approximately with tolerances
$\epsilon_{d,k}\searrow 0$ and $\epsilon_{p,k}\searrow 0$. This gives rise to
the analysis in Section 5.3.

###  5.3 Convergence of the outer iterations

For each $k\in\mathds{N}$, the stopping condition of Algorithm 2 occurs in a
finite number of iterations. Let $j_{k}$ denote the number of iterations
performed by Algorithm 2 at outer iteration $k$. To simplify the notation, let
$\bar{x}_{k}:=x_{k,j_{k}}$, $\bar{z}_{k}:=z_{k,j_{k}}$,
$\bar{s}_{k,1}:=s_{k,j_{k},1}$, $\bar{\Delta}_{k}:=\Delta_{k,j_{k}}$ and
$\bar{\nu}_{k}:=\nu_{k,j_{k}}$.

First, we present two assumptions that will be useful for our analysis.

######  Parameter Assumption 5.1.

$\liminf\Delta_{\min,k}=\bar{\Delta}>0$

To satisfy the second part of Parameter Assumption 5.1, we need to change
$\varphi$ in (27) to

(48) $\displaystyle\varphi(s;x_{k,j},B_{k,j})$
$\displaystyle:=\varphi_{f}(s;x_{k,j},B_{k,j})+\phi_{k}(x_{k,j})-\mu_{k}e^{T}X_{k,j}^{-1}s+\tfrac{1}{2}s^{T}\Theta_{k,j}s$
$\displaystyle=(f+\phi_{k})(x_{k,j})+(\nabla
f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e)^{T}s+\tfrac{1}{2}s^{T}(B_{k,j}+\Theta_{k,j})s,$

where $\Theta_{k,j}=\min(X_{k,j}^{-1}Z_{k,j},\kappa_{\textup{bar}}I)$ with the
$\min$ taken componentwise and $\kappa_{\textup{bar}}>0$.

Since the results of Section 5.2 involving $\varphi$ defined in (27) are all
based upon Step Assumption 5.1, those results continue to apply if Step
Assumption 5.1 holds for $\varphi$ defined in (48). We now show that that is
the case.

Fist, we observe that

$(f+\phi_{k})(x_{k,j}+s_{k,j})-\varphi(s_{k,j};x_{k,j},B_{k,j})=f(x_{k,j}+s_{k,j})-f(x_{k,j})-\nabla
f(x_{k,j})^{T}s_{k,j}+\\\
\phi_{k}(x_{k,j}+s_{k,j})-\phi_{k}(x_{k,j})-\mu_{k}s_{k,j}^{T}X_{k,j}^{-1}e-\tfrac{1}{2}s_{k,j}^{T}(B_{k,j}+\Theta_{k,j})s_{k,j}.$

Assume, as in Section 5.2, that $\|B_{k,j}\|\leq\tilde{M}_{k}$. Since
$\Theta_{k,j}$ is bounded by definition, we use (40) and (41) to conclude that
$f(x_{k,j}+s_{k,j})+\phi_{k}(x_{k,j}+s_{k,j})-\varphi(s_{k,j};x_{k,j},B_{k,j})=O(\|s_{k,j}\|^{2})$,
and, if $\psi(s;x_{k,j})=h(x_{k,j}+s)$, (35a) holds. The proof of Model
Assumption 5.2 is still valid when considering $B_{k,j}+\Theta_{k,j}$ instead
of $B_{k,j}+X_{k,j}^{-1}Z_{k,j}$, so that (35b) also holds. As a consequence,
Step Assumption 5.1 still holds with $\varphi$ defined in (48).

Now, our goal is to find a sequence $\\{w_{k}\\}\rightarrow 0$ such that
$w_{k}\in\nabla
f(\bar{x}_{k})-\bar{z}_{k}+\partial\psi(\bar{s}_{k,1};\bar{x}_{k})$. Under
some additional assumptions on $\psi$, this will allow us to establish that
Algorithm 1 generates iterates that satisfy asymptotically (9). We begin with
preliminary lemmas.

Assume that $s_{k,j,1}$ is not on the boundary of
$\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})$ and that
Model Assumption 5.1 holds. Then,

(49) $-\nu_{k,j}^{-1}s_{k,j,1}+\mu_{k}X_{k,j}^{-1}e\in\nabla
f(x_{k,j})+\partial\psi(s_{k,j,1};x_{k,j}).$

###### Proof.

The first step of the proximal gradient method $s_{k,j,1}$ satisfies (31).
According to (18), its first-order optimality conditions are

(50) $0\in
s_{k,j,1}+\nu_{k,j}\nabla\varphi_{\textup{cp}}(0;x_{k,j})+\partial(\nu_{k,j}\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})))(s_{k,j,1}).$

As $s_{k,j,1}$ is not on the boundary of
$\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})$, there exists
$r>0$ such that for all $s$ in an open ball of center $s_{k,j,1}$ and radius
$r$,
$\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))=0$.
Thus, the definition of the subdifferential guarantees that

$\partial(\nu_{k,j}\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})))(s_{k,j,1})=\\\
\nu_{k,j}\partial\psi(s_{k,j,1};x_{k,j})+\partial\chi(s_{k,j,1}\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})).$

We know that
$\partial\chi(s_{k,j,1}\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))=N_{\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})}(s_{k,j,1})$
using Section 5.1, and
$N_{\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})}(s_{k,j,1})=\\{0\\}$
because $s_{k,j,1}$ is not on the boundary of
$\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})$. Therefore,
(50) simplifies to

$0\in\nu_{k,j}^{-1}s_{k,j,1}+\nabla
f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e+\partial\psi(s_{k,j,1};x_{k,j})$

using Model Assumption 5.1 for $\nabla\varphi_{\textup{cp}}(0;x_{k,j})$.

Let Parameter Assumption 5.1 be satisfied. Then, for all $k\in\mathds{N}$,

$\|\bar{s}_{k,1}\|\leq\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}.$

###### Proof.

Since (33) holds, (45) leads to

$\epsilon_{d,k}\geq\bar{\nu}_{k}^{-1/2}\xi_{\textup{cp}}(\bar{\Delta}_{k};\bar{x}_{k},\bar{\nu}_{k})^{1/2}\geq\tfrac{1}{\sqrt{2}}\bar{\nu}_{k}^{-1}\|\bar{s}_{k,1}\|,$

which completes the proof.

The following assumption will be useful to establish that $\bar{s}_{k,1}$
converges to zero sufficiently fast to guarantee the convergence of the outer
iterations.

######  Parameter Assumption 5.2.

The sequences $\\{\epsilon_{d,k}\\}$ used in Algorithm 1 and
$\\{\kappa_{\textup{mdb},k}\\}$ from Section 5.2 satisfy

(51) $\kappa_{\textup{mdb},k}^{-1}\epsilon_{d,k}\rightarrow 0.$

To justify that Parameter Assumption 5.2 is reasonable, assume for simplicity
that $i$ is an index such that $(\bar{x}_{k})_{i}=\kappa_{\textup{mdb},k}$. We
have

$-e_{i}^{T}\bar{X}_{k}\bar{z}_{k}+\mu_{k}\leq\|\bar{X}_{k}\bar{z}_{k}-\mu_{k}e\|\leq\epsilon_{p,k},$

so that

$(\bar{x}_{k})_{i}(\bar{z}_{k})_{i}=e_{i}^{T}\bar{X}_{k}\bar{z}_{k}\geq\mu_{k}-\epsilon_{p,k},$

and, if $\epsilon_{p,k}<\mu_{k}$,

$\kappa_{\textup{mdb},k}^{-1}=\frac{1}{(\bar{x}_{k})_{i}}\leq\frac{(\bar{z}_{k})_{i}}{\mu_{k}-\epsilon_{p,k}}.$

We multiply the above inequality by $\epsilon_{p,k}$ to obtain

$\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}\leq\frac{(\bar{z}_{k})_{i}}{\mu_{k}\epsilon_{p,k}^{-1}-1}.$

If $\frac{\epsilon_{p,k}}{\mu_{k}}\rightarrow 0$, e.g.,
$\epsilon_{p,k}=\mu_{k}^{1+\gamma_{k}}$ with $0<\gamma_{k}<1$,
$\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}\rightarrow 0$. Then, the choice
$\epsilon_{d,k}=O(\epsilon_{p,k})$ guarantees that Parameter Assumption 5.2 is
satisfied.

Let Model Assumption 5.1, Parameter Assumption 5.1 and Parameter Assumption
5.2 be satisfied, and for all $j$, $(f+h)(x_{k,j})\geq(f+h)_{\textup{low},k}$.
Then, there exists $N\in\mathcal{N}_{\infty}$ such that for all $k\in N$,
$\bar{s}_{k,1}$ is not on the boundary of
$\bar{\Delta}_{k}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(\bar{x}_{k})$.

###### Proof.

For all $s\in\mathds{R}^{n}$, $\|s\|_{\infty}\leq\|s\|$, thus Section 5.3
leads to
$\|\bar{s}_{k,1}\|_{\infty}\leq\|\bar{s}_{k,1}\|\leq\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}$.
As a consequence, if $\sqrt{2}\epsilon_{d,k}\bar{\nu}_{k}<\Delta_{\min,k}$,
then $\bar{s}_{k,1}$ is not on the boundary of $\bar{\Delta}_{k}\mathds{B}$.
This is certainly true for $k$ sufficiently large, because $\bar{\Delta}>0$ in
Parameter Assumption 5.1 and $\epsilon_{d,k}\rightarrow 0$.

Now, we show that $\bar{s}_{k,1}$ is not on the boundary of
$\mathds{R}_{\delta_{k}}^{n}(\bar{x}_{k})$ if $k$ is large enough. First, we
point out that $\bar{\nu}_{k}\not\rightarrow\infty$, because

$\nu_{k,j}=\frac{1}{\|B_{k,j}\|+\|\Theta_{k,j}\|+\alpha^{-1}\Delta_{k,j}^{-1}}\leq\alpha\Delta_{k,j}.$

Then, we have

$\displaystyle
0\leq\frac{|\min_{i}(\bar{s}_{k,1})_{i}|}{\min_{i}(\bar{x}_{k})_{i}}\leq\frac{\|\bar{s}_{k,1}\|_{\infty}}{\min_{i}(\bar{x}_{k})_{i}}$
$\displaystyle\leq\frac{\|\bar{s}_{k,1}\|}{\min_{i}(\bar{x}_{k})_{i}}$
$\displaystyle\leq\frac{\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}}{\min_{i}(\bar{x}_{k})_{i}}$
using Section 5.3
$\displaystyle\leq\frac{\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}}{\kappa_{\textup{mdb},k}}$
$\displaystyle\text{using \lx@cref{creftypecap~refnum}{prop:min-xkji-kappa-
mdb}},$

and Parameter Assumption 5.2 indicates that
$\sqrt{2}\kappa_{\textup{mdb},k}^{-1}\bar{\nu}_{k}\epsilon_{d,k}\to 0$. As
$\delta_{k}\rightarrow\bar{\delta}<1$, the inequality

(52)
$\frac{\min_{i}(\bar{s}_{k,1})_{i}}{\min_{i}(\bar{x}_{k})_{i}}+1>\delta_{k}$

is satisfied if $k$ is large enough, and

(53) $\displaystyle\min_{i}(\bar{s}_{k,1}+\bar{x}_{k})_{i}$
$\displaystyle\geq\min_{i}(\bar{s}_{k,1})_{i}+\min_{i}(\bar{x}_{k})$ by
properties of the $\min$ $\displaystyle>\delta_{k}\min_{i}(\bar{x}_{k})$
$\displaystyle\text{ with~{}\eqref{eq:ineq-min-xskj}}.$

Therefore, $\bar{s}_{k,1}$ is not on the boundary of
$\mathds{R}_{\delta_{k}}(\bar{x}_{k})$ if $k$ is large enough. We conclude
that there exists $N\in\mathcal{N}_{\infty}$ such that for all $k\in N$,
$\bar{s}_{k,1}$ is not on the boundary of
$\bar{\Delta}_{k}\mathds{B}\cap\mathds{R}_{\delta_{k}}(\bar{x}_{k})$.

Let Model Assumption 5.1, Parameter Assumption 5.1 and Parameter Assumption
5.2 be satisfied, and for all $j$, $(f+h)(x_{k,j})\geq(f+h)_{\textup{low},k}$.
We define

(54)
$w_{k}:=-\bar{\nu}_{k}^{-1}\bar{s}_{k,1}+\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}.$

Then, there exists a subsequence $N\in\mathcal{N}_{\infty}$ such that for all
$k\in N$,

(55) $w_{k}\in\nabla
f(\bar{x}_{k})-\bar{z}_{k}+\partial\psi(\bar{s}_{k,1};\bar{x}_{k}),$

and

(56)
$\|w_{k}\|\leq\sqrt{2}\epsilon_{d,k}+\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}\rightarrow
0.$

###### Proof.

Section 5.3 indicates that there is a subsequence $N\in\mathcal{N}_{\infty}$
such that, for $k\in N$, $\bar{s}_{k,1}$ is not on the boundary of
$\bar{\Delta}_{k}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(\bar{x}_{k})$.
Thus, Section 5.3 holds, and

(57)
$-\bar{\nu}_{k}^{-1}\bar{s}_{k,1}+\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}\in\nabla
f(\bar{x}_{k})-\bar{z}_{k}+\partial\psi(\bar{s}_{k,1};\bar{x}_{k}).$

With $w_{k}$ defined in (54), we have (55). Now, for all
$i\in\\{1,\ldots,n\\}$, we use Section 5.2 to establish that

$(\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k})_{i}=(\mu_{k}-(\bar{x}_{k})_{i}(\bar{z}_{k})_{i})/(\bar{x}_{k})_{i}\leq(\mu_{k}-(\bar{x}_{k})_{i}(\bar{z}_{k})_{i})/\kappa_{\textup{mdb},k},$

and, by summing the square of the above inequality for all
$i\in\\{1,\ldots,n\\}$,

(58) $\displaystyle\|\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}\|^{2}$
$\displaystyle=\sum_{i=1}^{n}(\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k})_{i}^{2}$
$\displaystyle\leq\sum_{i=1}^{n}(\mu_{k}-(\bar{x}_{k})_{i}(\bar{z}_{k})_{i})^{2}/\kappa_{\textup{mdb},k}^{2}$
$\displaystyle=\|\mu_{k}e-\bar{X}_{k}\bar{z}_{k}\|^{2}/\kappa_{\textup{mdb},k}^{2}$
$\displaystyle\leq\kappa_{\textup{mdb},k}^{-2}\epsilon_{p,k}^{2}\text{
because~{}\eqref{eq:central-path-inner} holds},$

so that
$\|\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}\|\leq\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}$.
As $\epsilon_{p,k}\rightarrow 0$ and Parameter Assumption 5.2 holds, we deduce
that $\|\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}\|\rightarrow 0$.

Finally,

$\displaystyle\|w_{k}\|$
$\displaystyle=\|-\bar{\nu}_{k}^{-1}\bar{s}_{k,1}+\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}\|$
$\displaystyle\leq\bar{\nu}_{k}^{-1}\|\bar{s}_{k,1}\|+\|\mu_{k}\bar{X}_{k}^{-1}e-\bar{z}_{k}\|$
$\displaystyle\leq\sqrt{2}\epsilon_{d,k}+\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}\underset{k\rightarrow+\infty}{\longrightarrow}0,$

where we used Section 5.3 and (58) in the last inequality.

Now, we present two assumptions on $\psi$. The first will not be necessary for
the remaining results of this subsection, except as one of the justifications
for the second assumption. However, it will be used in Section 5.4.

######  Model Assumption 5.3.

$x\mapsto\psi(\cdot;x)$ is epi-continuous on $\mathds{R}_{+}^{n}$.

Model Assumption 5.3 holds if $\psi$ is continuous on
$\mathds{R}_{+}^{n}\times\mathds{R}_{+}^{n}$, but this condition is only
sufficient, not necessary [22, Exercise $7.40$]. Let us consider the case
where $\psi(\cdot;x)=s\mapsto h(x+s)$. Because $h$ is lsc, its epigraph is
closed, thus the sequence of functions $\\{h,h,\ldots\\}$ satisfies
$\\{h,h,\ldots\\}\overset{\textup{e}}{\longrightarrow}h$. Let
$\bar{x}\in\mathds{R}_{+}^{n}$ and $\\{x_{k}\\}\rightarrow\bar{x}$. [22,
Exercise $7.8d$] indicates that for $h_{k}=s\mapsto h(x_{k}+s)$ and
$\bar{h}=s\mapsto h(\bar{x}+s)$,
$h_{k}\overset{\textup{e}}{\longrightarrow}\bar{h}$. Since the latter is true
for all $\bar{x}\in\mathds{R}_{+}^{n}$, we conclude that Model Assumption 5.3
is satisfied.

######  Model Assumption 5.4.

For any sequences $\\{s_{k}\\}\rightarrow 0$ and
$\\{x_{k}\\}\rightarrow\bar{x}\geq 0$ such that $x_{k}+s_{k}>0$ for all
$k\in\mathds{N}$,

(59)
$\limsup_{k\rightarrow\infty}\partial\psi(s_{k};x_{k})\subset\partial\psi(0;\bar{x})=\partial
h(\bar{x}).$

We present some cases for which Model Assumption 5.4 holds.

  * •

When
$\mathop{\textup{g-lim}}_{k\rightarrow\infty}\partial\psi(\cdot;x_{k})=\partial\psi(\cdot;\bar{x})$.
Attouch’s theorem [3] (also written in [22, Theorem $12.35$]) indicates that
this condition is satisfied when $\psi(\cdot;x_{k})$ and $\psi(\cdot;\bar{x})$
are proper, lsc, convex functions with
$\psi(\cdot;x_{k})\overset{\textup{e}}{\longrightarrow}\psi(\cdot;\bar{x})$
(i.e., Model Assumption 5.3 holds). An extension to non-convex functions under
some more sophisticated assumptions is established by Poliquin [21].

  * •

When $\psi(s;x)=h(x+s)$ and $h(x_{k}+s_{k})\rightarrow h(\bar{x})$ (e.g.,
$h=\|\cdot\|_{1}$), using Section 2 applied to $\\{x_{k}+s_{k}\\}\to\bar{x}$.

  * •

When $\psi(s;x)=h(x+s)$ but $h$ is not continuous, we may still be able to
show that Model Assumption 5.4 holds. For example, with $h=\|\cdot\|_{0}$,
[16, Theorem $1$] shows that $\partial h(x)=\\{v\mid v_{i}=0\text{ if
}x_{i}\neq 0\\}$. Thus,
$\partial\psi(s_{k};x_{k})=\partial(\|\cdot\|_{0})(x_{k}+s_{k})=\\{0\\}\subset\partial
h(\bar{x})$.

Under the assumptions of Section 5.3 and Model Assumption 5.4, let $\bar{x}$
and $\bar{z}$ be limit points of $\\{\bar{x}_{k}\\}$ and $\\{\bar{z}_{k}\\}$,
respectively. Then,

(60) $0\in\nabla f(\bar{x})-\bar{z}+\partial
h(\bar{x})\quad\text{and}\quad\bar{X}\bar{Z}e=0.$

In this case, when the CQ is satisfied at $\bar{x}$, $\bar{x}$ is first-order
stationary for (1).

###### Proof.

In Section 5.3, $N$ can be chosen such that $\bar{x}_{k}\underset{k\in
N}{\longrightarrow}\bar{x}$ and $\bar{z}_{k}\underset{k\in
N}{\longrightarrow}\bar{z}$. We apply Model Assumption 5.4 to (55) and (56) to
deduce

(61) $-\nabla f(\bar{x})+\bar{z}\in\partial\psi(0;\bar{x})=\partial
h(\bar{x}),$

which indicates that $0\in
P^{\mathcal{L}}(\bar{\Delta};\bar{x},\bar{z},\bar{\nu})$. The condition (34)
implies that $\bar{X}\bar{Z}e=0$. When the CQ is satisfied, we can use Section
3 to conlude that $\bar{x}$ is first-order stationary for (1).

In Section 5.3, $w_{k}$ evokes of the concept of
$(\epsilon_{p},\epsilon_{d})$-KKT optimality for interior-point methods
introduced in [11, Definition $2.1$]. We slightly modify this concept in the
following definition.

######  Definition 5.1.

Let $\epsilon_{p},\epsilon_{d}\geq 0$, $x\geq 0$, and $\Delta,\nu>0$. $x$ is
said to be $(\epsilon_{p},\epsilon_{d})$-KKT optimal if there exist $z\geq 0$
and $v\in\partial h(x)$ such that

(62) $\|\nabla f(x)-z+v\|\leq\epsilon_{d},$

and, for all $i\in\\{1,\ldots,n\\}$,

(63) $x_{i}z_{i}\leq\epsilon_{p}.$

The main modification to the original formulation in [11, Definition $2.1$] is
that we require $x_{i}z_{i}\leq\epsilon_{p}$ instead of
$\min(x_{i},z_{i})\leq\epsilon_{p}$ for all $i$, but this is linked to our
different choices of stopping condition for the complementary slackness. The
first part of the definition (62) is similar to the $\epsilon$-stationarity
[10, Definition $4.5$] for more general problems.

Model Assumption 5.1 does not necessarily guarantee that
$\partial\psi(s_{k,j,1};x_{k,j})=\partial h(x_{k,j}+s_{k,j,1})$. Thus, for
$k\in N$ where $N$ is a subsequence introduced in Section 5.3, we cannot use
Section 5.3 to measure the $(\epsilon_{p},\epsilon_{d})$-KKT optimality of
$\bar{x}_{k}$. Let

(64) $\epsilon_{h,k}=\mathop{\textup{dist}}(w_{k}-\nabla
f(\bar{x}_{k})+\bar{z}_{k},\,\partial h(\bar{x}_{k}+\bar{s}_{k,1})).$

As Section 5.3 indicates that
$w_{k}-f(\bar{x}_{k})+\bar{z}_{k}\in\partial\psi(\bar{s}_{k,1};\bar{x}_{k})$,
we can obtain a measure of $(\epsilon_{p},\epsilon_{d})$-KKT optimality which
depends on $\epsilon_{h,k}$. When all the elements of
$\partial\psi(\bar{s}_{k,1};\bar{x}_{k})$ are close to an element of $\partial
h(\bar{x}_{k}+\bar{s}_{k,1})$, we expect $\epsilon_{h,k}$ to be small. In
particular, if $\partial\psi(\bar{s}_{k,1};\bar{x}_{k})\subseteq\partial
h(\bar{x}_{k}+\bar{s}_{k,1})$, $\epsilon_{h,k}=0$. Let the assumptions of
Section 5.3 be satisfied, and $\epsilon_{h,k}$ be defined in (64). Then, there
exists $N\in\mathcal{N}_{\infty}$ such that for all $k\in N$,
$\bar{x}_{k}+\bar{s}_{k,1}$ is
$(\bar{\epsilon}_{p,k},\bar{\epsilon}_{d,k})$-KKT optimal with constants

$\displaystyle\bar{\epsilon}_{p,k}$
$\displaystyle=\epsilon_{p,k}+\sqrt{n}\mu_{k}+\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}\|\bar{z}_{k}\|$
$\displaystyle\bar{\epsilon}_{d,k}$
$\displaystyle=\epsilon_{h,k}+\sqrt{2}\epsilon_{d,k}(1+\bar{\nu}_{k}L_{f})+\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}.$

######  Proof 5.2.

Section 5.3 guarantees that, for all $k$ in a subsequence
$N\in\mathcal{N}_{\infty}$, (55) holds, i.e.

$w_{k}-\nabla
f(\bar{x}_{k})+\bar{z}_{k}\in\partial\psi(\bar{s}_{k,1};\bar{x}_{k}),$

where $w_{k}$ is defined in (54). As $\partial h(\bar{x}_{k}+\bar{s}_{k,1})$
is closed, we can choose $v_{k}\in\partial h(\bar{x}_{k}+\bar{s}_{k,1})$ such
that, for $y_{k}:=v_{k}-w_{k}+\nabla f(\bar{x}_{k})-\bar{z}_{k}$, we have
$\|y_{k}\|=\epsilon_{h,k}$. Then,

$(v_{k}-w_{k}+\nabla f(\bar{x}_{k})-\bar{z}_{k})+w_{k}\in\nabla
f(\bar{x}_{k})-\bar{z}_{k}+\partial h(\bar{x}_{k}+\bar{s}_{k,1}),$

which we may rewrite as

$y_{k}+w_{k}+\nabla f(\bar{x}_{k}+\bar{s}_{k,1})-\nabla
f(\bar{x}_{k})\in\nabla f(\bar{x}_{k}+\bar{s}_{k,1})-\bar{z}_{k}+\partial
h(\bar{x}_{k}+\bar{s}_{k,1}),$

and

$\|y_{k}+w_{k}+\nabla f(\bar{x}_{k}+\bar{s}_{k,1})-\nabla
f(\bar{x}_{k})\|\leq\epsilon_{h,k}+\|w_{k}\|+L_{f}\|\bar{s}_{k,1}\|.$

Section 5.3 implies that
$L_{f}\|\bar{s}_{k,1}\|\leq\sqrt{2}L_{f}\bar{\nu}_{k}\epsilon_{d,k}$. We
combine the latter inequality with (56) in Section 5.3 to obtain

(65) $\|y_{k}+w_{k}+\nabla f(\bar{x}_{k}+\bar{s}_{k,1})-\nabla
f(\bar{x}_{k})\|\leq\epsilon_{h,k}+\sqrt{2}\epsilon_{d,k}(1+\bar{\nu}_{k}L_{f})+\kappa_{\textup{mdb},k}^{-1}\epsilon_{p,k}.$

Now, for all $i\in\\{1,\ldots,n\\}$,

(66) $\displaystyle(\bar{x}_{k}+\bar{s}_{k,1})_{i}(\bar{z}_{k})_{i}$
$\displaystyle\leq\|(\bar{X}_{k}+\bar{S}_{k,1})\bar{z}_{k}\|$
$\displaystyle\leq\|\bar{X}_{k}\bar{z}_{k}-\mu_{k}e\|+\|\mu_{k}e\|+\|\bar{S}_{k,1}\bar{z}_{k}\|$
$\displaystyle\leq\epsilon_{p,k}+\sqrt{n}\mu_{k}+\|\bar{s}_{k,1}\|\|\bar{z}_{k}\|$
$\displaystyle\leq\epsilon_{p,k}+\sqrt{n}\mu_{k}+\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}\|\bar{z}_{k}\|.$

We use (65), (66) and Definition 5.1 to conclude.

If $\bar{z}_{k}$ is bounded,
$\bar{\epsilon}_{p,k}\underset{N}{\longrightarrow}0$ in Equation 64. If
$\epsilon_{h,k}\rightarrow 0$, we also have
$\bar{\epsilon}_{d,k}\underset{N}{\longrightarrow}0$.

###  5.4 Convergence with a new criticality measure

Now, instead of using
$\nu_{k,j}^{-1/2}\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}$
(involving $\varphi_{\textup{cp}}(\cdot;x_{k,j})$) for the criticality measure
of Algorithm 2, we would like to use a measure based upon
$\varphi^{\mathcal{L}}(\cdot;x_{k,j},z_{k,j})$ defined in (10a). The reason
behind this choice is inspired from the criticality measure $\|\nabla
f(x_{k,j})-z_{k,j}\|$ used in primal-dual trust region algorithms in the
smooth case, instead of $\|\nabla f(x_{k,j})-\mu_{k}X_{k,j}^{-1}e\|$ used in
primal algorithms, see for example [9, Algorithm $13.6.2$]. We may expect that
this choice results in fewer iterations of Algorithm 2 when $x_{k,j}$ and
$z_{k,j}$ are close to a solution of (1), because (we express this idea with
smooth notations for now), if $j_{k}$ is the index for which the stopping
criteria of Algorithm 2 are met, $\|\nabla
f(\bar{x}_{k})-\bar{z}_{k}\|=\|\nabla f(x_{k+1,0})-z_{k+1,0}\|$, whereas
$\|\nabla f(\bar{x}_{k})-\mu_{k}\bar{X}_{k}^{-1}e\|\neq\|\nabla
f(x_{k+1,0})-\mu_{k+1}X_{k+1,0}^{-1}e\|$. However, to change the stopping
criterion, we need the following convexity assumption.

######  Model Assumption 5.5.

For a sequence $\\{x_{k,j}\\}_{j}$ generated by Algorithm 2 at iteration $k$,
$\psi(\cdot;x_{k,j})$ is convex for all $j$.

If $\psi(s;x)=h(x+s)$ and $h$ is convex, Model Assumption 5.5 holds.

In this section, we define

(67)
$s_{k,j}^{\mathcal{L}}\in\mathop{\textup{argmin}}_{s}m^{\mathcal{L}}(s;x_{k,j},z_{k,j},\nu_{k,j})+\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})),$

where $m^{\mathcal{L}}(s;x_{k,j},z_{k,j}\nu_{k,j})$ is defined in (11), and

(68)
$\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,j};x_{k,j},z_{k,j},\nu_{k,j})=(f+\phi_{k}+h)(x_{k,j})-(\varphi^{\mathcal{L}}(s_{k,j}^{\mathcal{L}};x_{k,j},z_{k,j})-\psi(s_{k,j}^{\mathcal{L}};x_{k,j})),$

where $\varphi^{\mathcal{L}}(\cdot;x_{k,j},z_{k,j})$ is defined in (10a). We
point out that $\xi_{\delta_{k}}^{\mathcal{L}}$ and $\xi^{\mathcal{L}}$
defined in (13) are almost identical, the latter being computed by replacing
$\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))$ by
$\chi(s\mid\Delta_{k,j}\mathds{B})+\chi(x_{k,j}+s\mid\mathds{R}_{+}^{n})$ in
(67). Model Assumption 5.5 will be useful in Section 5.4, which is crucial for
our analysis of Algorithm 3 because it establishes the convergence of the
inner iterations with $\xi_{\delta_{k}}^{\mathcal{L}}$.

Algorithm 3 resembles Algorithm 1, except for the stopping criterion (33). Our
ultimate goal in this subsection is to show that replacing (33) in Algorithm 1
by (69) in Algorithm 3 maintains similar convergence properties to those of
Section 5.3, but, as illustrated in Section 6, performs better in practice.

Algorithm 3 Nonsmooth interior-point method (outer iteration) with stopping
criteria based upon $\xi_{\delta_{k}}^{\mathcal{L}}$ in (68).

1:Choose $\epsilon>0$, sequences $\\{\mu_{k}\\}\searrow 0$,
$\\{\epsilon_{d,k}\\}\searrow 0$, $\\{\epsilon_{p,k}\\}\searrow 0$, and
$\\{\delta_{k}\\}\rightarrow\bar{\delta}\in[0,1)$ with $\delta_{k}\in(0,1)$
for all $k$.

2:Choose $x_{0,0}\in\mathds{R}^{n}_{++}$ where $h$ is finite.

3:for $k=0,1,\ldots$ do

4: Compute an approximate solution $x_{k}:=x_{k,j}$ to (22) and
$z_{k}:=z_{k,j}$ in the sense that  (69)
$\nu_{k,j}^{-1/2}\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,j};x_{k,j},\nu_{k,j})^{1/2}\leq\epsilon_{d,k},$
and (34) holds, that we recall in the following inequality for conveniency
$\|\mu_{k}e-X_{k,j}z_{k,j}\|\leq\epsilon_{p,k}.$

5: Set $x_{k+1,0}:=x_{k}$.

6:end for

The following result shows that the inner iteration terminates finitely.

Under the assumptions of Section 5.2, Model Assumption 5.3 and Model
Assumption 5.5, if $\\{x_{k,j}\\}_{j}$ possesses a limit point $x_{k}^{*}$,
then

(70)
$\liminf_{j\rightarrow+\infty}\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,j};x_{k,j},z_{k,j},\nu_{k,j})=0,$

and

(71) $\liminf_{j\rightarrow+\infty}\|s_{k,j}^{\mathcal{L}}\|=0,$

where $s_{k,j}^{\mathcal{L}}$ is defined in (67).

######  Proof 5.3.

As $\Delta_{k,j}\in[\Delta_{\min,k},\Delta_{\max}]$ and
$\nu_{k,j}\in[\nu_{\min,k},1]$, there exists an infinite subsequence $N$ such
that $\Delta_{k,j}\underset{j\in N}{\longrightarrow}\Delta_{k}^{*}$,
$\nu_{k,j}\underset{j\in N}{\longrightarrow}\nu_{k}^{*}$,
$x_{k,j}\underset{j\in N}{\longrightarrow}x_{k}^{*}$, and with Section 5.2
$z_{k,j}\underset{j\in N}{\longrightarrow}z_{k}^{*}$ with
$X_{k}^{*}Z_{k}^{*}e=\mu_{k}e$. By continuity of the $\min$,
$\mathds{R}_{\delta_{k}}^{n}(x_{k,j})\underset{j\in
N}{\longrightarrow}\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*})$. The sets
$\Delta_{k,j}\mathds{B}$ and $\mathds{R}_{\delta_{k}}^{n}(x_{k,j})$ are convex
(using Section 5.1 for the latter). Since $\Delta_{k}^{*}\mathds{B}$ and
$\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*})$ are convex and cannot be separated,
we use [22, Theorem $4.33$] to conclude that

$\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})\underset{j\in
N}{\longrightarrow}\Delta_{k}^{*}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*}).$

With [22, Theorem $7.4f$], we deduce

$\mathop{\textup{e-lim}}_{j\in
N}\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))=\chi(\cdot\mid\Delta_{k}^{*}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*})).$

Thanks to Section 5.2 and the smoothness of $\varphi_{\textup{cp}}(\cdot;x)$
and $\varphi^{\mathcal{L}}(\cdot;x,z)$, we also have

(72) $\mathop{\textup{e-lim}}_{j\in
N}\varphi_{\textup{cp}}(\cdot;x_{k,j})=\mathop{\textup{e-lim}}_{j\in
N}\varphi^{\mathcal{L}}(\cdot;x_{k,j},z_{k,j})=\varphi^{\mathcal{L}}(\cdot;x_{k}^{*},z_{k}^{*}).$

The functions $\varphi_{\textup{cp}}(\cdot;x_{k,j})$,
$\varphi^{\mathcal{L}}(\cdot;x_{k,j},z_{k,j})$ and
$\varphi^{\mathcal{L}}(\cdot;x_{k}^{*},z_{k}^{*})$ are all convex because they
are linear. Model Assumption 5.3 implies that $\mathop{\textup{e-lim}}_{j\in
N}\psi(\cdot;x_{k,j})=\psi(\cdot;x_{k}^{*})$. Model Assumption 5.5 and [22,
Theorem $7.46$] lead to

(73) $\mathop{\textup{e-lim}}_{j\in
N}\psi(\cdot;x_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))=\psi(\cdot;x_{k}^{*})+\chi(\cdot\mid\Delta_{k}^{*}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*})),$

and the above functions are all convex. We deduce from (72), (73) and [22,
Theorem $7.46$] that

$\mathop{\textup{e-lim}}_{j\in
N}m_{\textup{cp}}(\cdot;x_{k,j},\nu_{k,j})=\mathop{\textup{e-lim}}_{j\in
N}m^{\mathcal{L}}(\cdot;x_{k,j},z_{k,j},\nu_{k,j})=m^{\mathcal{L}}(\cdot;x_{k}^{*},z_{k}^{*},\nu_{k}^{*}),$

where $m^{\mathcal{L}}$ is defined in (11). The sequences
$m_{\textup{cp}}(\cdot;x_{k,j},\nu_{k,j})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))$
and
$m^{\mathcal{L}}(\cdot;x_{k}^{*},z_{k}^{*},\nu_{k}^{*})+\chi(\cdot\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))$
are level-bounded because of the indicators. As in (31), we have

$\mathop{\textup{prox}}_{\nu_{k}^{*}\psi(\cdot;x_{k}^{*})+\chi(\cdot\mid\Delta_{k}^{*}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*}))}(-\nu_{k}^{*}\nabla\varphi_{\textup{cp}}(0;x_{k}^{*},\nu_{k}^{*}))=\\\
\arg\min_{s}m^{\mathcal{L}}(s;x_{k}^{*},z_{k}^{*},\nu_{k}^{*})+\chi(s\mid\Delta_{k}^{*}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k}^{*})),$

and the above problem is single valued because of [22, Theorem $2.26a$]. Let
$s^{*}$ denote its only solution. Section 2, and specifically (5), implies
that the sequences

$s_{k,j,1}\in\arg\min_{s}m_{\textup{cp}}(s;x_{k,j},\nu_{k,j})+\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})),$

and

$s_{k,j}^{\mathcal{L}}\in\arg\min_{s}m^{\mathcal{L}}(s;x_{k,j},z_{k,j},\nu_{k,j})+\chi(s\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))$

have the same limit $s^{*}$. We have shown in Equation 45 that
$s_{k,j,1}\underset{j\rightarrow+\infty}{\longrightarrow}0$. Thus, $s^{*}=0$.
Finally, we have

$\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})=(f+h+\phi_{k})(x_{k,j})-m_{\textup{cp}}(s_{k,j,1};x_{k,j},\nu_{k,j})+\tfrac{1}{2}\nu_{k,j}^{-1}\|s_{k,j,1}\|^{2}.$

As
$\xi_{\textup{cp}}(\Delta_{k,j};x_{k,j},\nu_{k,j})\underset{j\rightarrow+\infty}{\longrightarrow}0$
by Model Assumption 5.2, and
$\|s_{k,j,1}\|\underset{j\rightarrow+\infty}{\longrightarrow}0$, we deduce
that

$m_{\textup{cp}}(s_{k,j,1};x_{k,j},\nu_{k,j})\underset{j\in
N}{\longrightarrow}(f+h+\phi_{k})(x_{k}^{*}).$

Using (4) in Section 2, we have

(74) $\lim_{j\in N}m_{\textup{cp}}(s_{k,j,1};x_{k,j},\nu_{k,j})=\lim_{j\in
N}m^{\mathcal{L}}(s_{k,j}^{\mathcal{L}};x_{k,j},z_{k,j},\nu_{k,j})=(f+h+\phi_{k})(x_{k}^{*}).$

The expression of $\xi_{\delta_{k}}^{\mathcal{L}}$ in (68) can also be written
as

$\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,j};x_{k,j},z_{k,j},\nu_{k,j})=(f+h+\phi_{k})(x_{k,j})-m^{\mathcal{L}}(s_{k,j}^{\mathcal{L}};x_{k,j},z_{k,j},\nu_{k,j})+\tfrac{1}{2}\nu_{k,j}^{-1}\|s_{k,j}^{\mathcal{L}}\|^{2}.$

By injecting the limit of $s_{k,j}^{\mathcal{L}}$ and (74) in the above
equation, we obtain (70).

From this point on, $j_{k}$ denotes the number of iterations performed by
Algorithm 2 at iteration $k$ with the inner stopping criteria from Algorithm
3, and we use again the notation $\bar{x}_{k}=x_{k,j_{k}}$,
$\bar{z}_{k}=z_{k,j_{k}}$, $\bar{s}_{k,1}=s_{k,j_{k},1}$,
$\bar{\Delta}_{k}=\Delta_{k,j_{k}}$, $\bar{\nu}_{k}=\nu_{k,j_{k}}$, with the
addition of $\bar{s}_{k}^{\mathcal{L}}:=s_{k,j_{k}}^{\mathcal{L}}$. The
following three lemmas are analogous to Section 5.3, Section 5.3 and Section
5.3.

Assume that $s_{k,j}^{\mathcal{L}}$ is not on the boundary of
$\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})$ and that
Model Assumption 5.1 holds. Then,

(75) $-\nu_{k,j}^{-1}s_{k,j}^{\mathcal{L}}\in\nabla
f(x_{k,j})-z_{k,j}+\partial\psi(s_{k,j}^{\mathcal{L}};x_{k,j}).$

######  Proof 5.4.

The first-order stationarity condition of (67) is

$-\nu_{k,j}^{-1}s_{k,j}^{\mathcal{L}}\in\nabla
f(x_{k,j})-z_{k,j}+\partial\psi(s_{k,j}^{\mathcal{L}};\bar{x}_{k})+\partial\chi(s_{k,j}^{\mathcal{L}}\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j})).$

The same analysis as in the proof of Section 5.3 establishes that

$\partial\chi(s_{k,j}^{\mathcal{L}}\mid\Delta_{k,j}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(x_{k,j}))=\\{0\\},$

so that (75) holds.

Let Parameter Assumption 5.1 be satisfied. Then, for all $k\in\mathds{N}$,

$\|\bar{s}_{k}^{\mathcal{L}}\|\leq\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}.$

######  Proof 5.5.

The bound

(76)
$\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,j};x_{k,j},z_{k,j},\nu_{k,j})\geq\tfrac{1}{2}\nu_{k,j}^{-1}\|s_{k,j}^{\mathcal{L}}\|^{2}$

holds because

$\displaystyle m^{\mathcal{L}}(0;x_{k,j},z_{k,j},\nu_{k,j})$
$\displaystyle=(f+\phi_{k}+h)(x_{k,j})$ $\displaystyle\geq
m^{\mathcal{L}}(s_{k,j}^{\mathcal{L}};x_{k,j},z_{k,j},\nu_{k,j})$
$\displaystyle=(f+\phi_{k}+h)(x_{k,j})+(\nabla
f(x_{k,j})+z_{k,j})^{T}s_{k,j}^{\mathcal{L}}+\tfrac{1}{2}\nu_{k,j}^{-1}\|s_{k,j}^{\mathcal{L}}\|^{2}.$

The stopping criterion (69) and (76) lead to

$\tfrac{1}{\sqrt{2}}\bar{\nu}_{k}^{-1}\|\bar{s}_{k}^{\mathcal{L}}\|\leq\bar{\nu}_{k}^{-1/2}\xi_{\delta_{k}}^{\mathcal{L}}(\bar{\Delta}_{k};\bar{x}_{k},\bar{z}_{k},\bar{\nu}_{k})^{1/2}\leq\epsilon_{d,k},$

which completes the proof.

Let Model Assumption 5.1, Parameter Assumption 5.1 and Parameter Assumption
5.2 be satisfied, and for all $j$, $(f+h)(x_{k,j})\geq(f+h)_{\textup{low},k}$.
Then, there exists $N\in\mathcal{N}_{\infty}$ such that for all $k\in N$,
$\bar{s}_{k}^{\mathcal{L}}$ is not on the boundary of
$\bar{\Delta}_{k}\mathds{B}\cap\mathds{R}_{\delta_{k}}^{n}(\bar{x}_{k})$.

######  Proof 5.6.

Since for any $s\in\mathds{R}^{n}$, $\|s\|_{\infty}\leq\|s\|$, Section 5.4
leads to
$\|\bar{s}_{k}^{\mathcal{L}}\|_{\infty}\leq\|\bar{s}_{k}^{\mathcal{L}}\|\leq\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}$.
Thus, if $\sqrt{2}\epsilon_{d,k}\bar{\nu}_{k}<\Delta_{\min,k}$,
$\bar{s}_{k}^{\mathcal{L}}$ is not on the boundary of
$\bar{\Delta}_{k}\mathds{B}$. As $\bar{\Delta}>0$ in Parameter Assumption 5.1
and $\epsilon_{d,k}\rightarrow 0$, this is true if $k$ is sufficiently large.
The rest of the proof is identical to that of Section 5.3.

Now, we can establish results similar to Section 5.3, Section 5.3, and
Equation 64 for Algorithm 3.

Let Model Assumption 5.1, Parameter Assumption 5.1 and Parameter Assumption
5.2 be satisfied, and for all $j$, $(f+h)(x_{k,j})\geq(f+h)_{\textup{low},k}$.
Then, there exists a subsequence $N\in\mathcal{N}_{\infty}$ such that for all
$k\in N$,

(77) $-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}\in\nabla
f(\bar{x}_{k})-\bar{z}_{k}+\partial\psi(\bar{s}_{k}^{\mathcal{L}};\bar{x}_{k}),$

with

(78)
$\nu_{k,j}^{-1}\|\bar{s}_{k}^{\mathcal{L}}\|\leq\sqrt{2}\epsilon_{d,k}\rightarrow
0.$

######  Proof 5.7.

Section 5.4 and Section 5.4 lead to (77). Section 5.4 shows that
$\nu_{k,j}^{-1}\|s_{k,j}^{\mathcal{L}}\|\leq\sqrt{2}\epsilon_{d,k}$, and, as
$\epsilon_{d,k}\rightarrow 0$, (78) is satisfied.

Under the assumptions of Section 5.4 and Model Assumption 5.4, let $\bar{x}$
and $\bar{z}$ be limit points of $\\{\bar{x}_{k}\\}$ and $\\{\bar{z}_{k}\\}$,
respectively. Then,

(79) $0\in\nabla f(\bar{x})-\bar{z}+\partial
h(\bar{x})\quad\text{and}\quad\bar{X}\bar{Z}e=0.$

In this case, when the CQ is satisfied, $\bar{x}$ is first-order stationary
for (1).

######  Proof 5.8.

In Section 5.4, $N$ can be chosen such that $\bar{x}_{k}\underset{k\in
N}{\longrightarrow}\bar{x}$ and $\bar{z}_{k}\underset{k\in
N}{\longrightarrow}\bar{z}$. We apply Model Assumption 5.4 to (77) and (78) to
obtain

$-\nabla f(\bar{x})+\bar{z}\in\partial\psi(0;\bar{x})=\partial h(\bar{x}),$

and we conclude as in the proof of Section 5.3.

Finally, we show a result similar to Equation 64 for the
$(\epsilon_{p},\epsilon_{d})$-KKT optimality. We use again a subsequence $N$
as in section 5.4, $v_{k}\in\partial
h(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})$, and we define

(80)
$\epsilon_{h,k}^{\mathcal{L}}=\mathop{\textup{dist}}(-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}-\nabla
f(\bar{x}_{k})+\bar{z}_{k},\,\partial h(\bar{x}_{k}+\bar{s}_{k,1})).$

Let the assumptions of Section 5.4 be satisfied, and
$\epsilon_{h,k}^{\mathcal{L}}$ be defined in (80). Then, there exists a
subsequence $N\in\mathcal{N}_{\infty}$ such that for all $k\in N$,
$\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}}$ is
$(\bar{\epsilon}_{p,k}^{\mathcal{L}},\bar{\epsilon}_{d,k}^{\mathcal{L}})$-KKT
optimal with constants

$\displaystyle\bar{\epsilon}_{p,k}^{\mathcal{L}}$
$\displaystyle=\epsilon_{p,k}+\sqrt{n}\mu_{k}+\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}\|\bar{z}_{k}\|$
$\displaystyle\bar{\epsilon}_{d,k}^{\mathcal{L}}$
$\displaystyle=\epsilon_{h,k}^{\mathcal{L}}+\sqrt{2}\epsilon_{d,k}(1+L_{f}\bar{\nu}_{k}).$

######  Proof 5.9.

Section 5.4 guarantees that there exists an infinite subsequence $N$ such that
for all $k\in N$,

$-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}-\nabla
f(\bar{x}_{k})+\bar{z}_{k}\in\partial\psi(\bar{s}_{k}^{\mathcal{L}};\bar{x}_{k}),$

Since $\partial h(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})$ is closed and
nonempty, we can choose $v_{k}\in\partial
h(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})$ such that, for
$y_{k}^{\mathcal{L}}:=-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}-\nabla
f(\bar{x}_{k})+\bar{z}_{k}$, we have
$\|y_{k}^{\mathcal{L}}\|=\epsilon_{h,k}^{\mathcal{L}}$. Now,

$(v_{k}-(-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}-\nabla
f(\bar{x}_{k})+\bar{z}_{k}))-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}\in\nabla
f(\bar{x}_{k})-\bar{z}_{k}+\partial h(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}}),$

which can also be written as

$y_{k}^{\mathcal{L}}-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}+\nabla
f(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})-\nabla f(\bar{x}_{k})\in\nabla
f(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})-\bar{z}_{k}+\partial
h(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}}).$

The triangle inequality combined and the Lipschitz constant $L_{f}$ of $\nabla
f$ leads to

$\|y_{k}^{\mathcal{L}}-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}+\nabla
f(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})-\nabla
f(\bar{x}_{k})\|\leq\epsilon_{h,k}^{\mathcal{L}}+\bar{\nu}_{k}^{-1}\|\bar{s}_{k}^{\mathcal{L}}\|+L_{f}\|\bar{s}_{k}^{\mathcal{L}}\|.$

Section 5.4 then implies that

$\|y_{k}^{\mathcal{L}}-\bar{\nu}_{k}^{-1}\bar{s}_{k}^{\mathcal{L}}+\nabla
f(\bar{x}_{k}+\bar{s}_{k}^{\mathcal{L}})-\nabla
f(\bar{x}_{k})\|\leq\epsilon_{h,k}^{\mathcal{L}}+\sqrt{2}\epsilon_{d,k}(1+L_{f}\bar{\nu}_{k}).$

Finally, the inequalities of (66) still hold when replacing $\bar{s}_{k,1}$ by
$\bar{s}_{k}^{\mathcal{L}}$: for all $i\in\\{1,\ldots,n\\}$,
$(\bar{x}_{k})_{i}(\bar{z}_{k})_{i}\leq\epsilon_{p,k}+\sqrt{n}\mu_{k}+\sqrt{2}\bar{\nu}_{k}\epsilon_{d,k}\|\bar{z}_{k}\|$.

When $\bar{x}_{k}\underset{k\in
N}{\longrightarrow}\bar{x}\in\mathds{R}_{+}^{n}$ and
$\bar{z}_{k}\underset{k\in N}{\longrightarrow}\bar{z}\in\mathds{R}_{+}^{n}$,
$\\{\bar{z}_{k}\\}_{N}$ is bounded, thus
$\bar{\epsilon}_{p,k}^{\mathcal{L}}\underset{N}{\longrightarrow}0$ in Section
5.4. If $\epsilon_{h,k}^{\mathcal{L}}\rightarrow 0$, we also have
$\bar{\epsilon}_{d,k}^{\mathcal{L}}\underset{N}{\longrightarrow}0$.

##  6 Implementation and numerical experiments

All solvers tested are available from RegularizedOptimization.jl. We define

$\epsilon_{d,k}:=\epsilon_{k}+\epsilon_{r,i}\nu_{k,0}^{-1/2}\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,0};x_{k,0},z_{k,0},\nu_{k,0})^{1/2}\quad\text{and}\quad\epsilon_{p,k}:=\epsilon_{k},$

where $\epsilon_{r,i}\geq 0$ is a predefined relative tolerance for the inner
iterations. Algorithm 2 terminates when

$\displaystyle\nu_{k,j}^{-1/2}\xi_{\delta_{k}}^{\mathcal{L}}(\Delta_{k,j};x_{k,j},z_{k,j},\nu_{k,j_{k}})^{1/2}$
$\displaystyle<\epsilon_{d,k}$ $\displaystyle\|X_{k,j}z_{k,j}-\mu_{k}e\|$
$\displaystyle<\epsilon_{p,k}.$

We use the constant $\kappa_{\textup{bar}}=10^{6}$ for $\Theta_{k,j}$ in (48)
and $\epsilon_{r,i}=10^{-1}$. We set $\epsilon_{0}=\mu_{0}=1$,
$\mu_{k+1}=\mu_{k}/10$, $\epsilon_{k}=\mu_{k}^{1.01}$, and
$\Delta_{k,0}=1000\mu_{k}$, similarly as Conn et al. [8] did in the smooth
case for their interior-point trust-region algorithm. In addition, we go to
iteration $k+1$ and set $x_{k+1,0}=x_{k,j}$ if Algorithm 2 performs more than
$j=200$ iterations. Any iteration $k$ where $x_{k,0}$ is not first-order
stationary for (1) has $\epsilon_{d,k}>\epsilon_{p,k}=\epsilon_{k}$, thus the
remarks below Parameter Assumption 5.2 are not sufficient to prove that this
assumption is always satisfied, however, we observe satisfying performance
with these parameters. Although it is possible to use
$\epsilon_{d,k}=\epsilon_{k}$, it would require more inner iterations.

To declare convergence of Algorithm 3, we use the following criteria

(81a) $\displaystyle\mu_{k}$ $\displaystyle<\epsilon$ (81b)
$\displaystyle\|\bar{X}_{k}\bar{z}_{k}-\mu_{k}e\|$ $\displaystyle<\epsilon$
(81c)
$\displaystyle\bar{\nu}_{k}^{-1/2}\xi_{\delta_{k}}^{\mathcal{L}}(\bar{\Delta}_{k};\bar{x}_{k},\bar{z}_{k},\bar{\nu}_{k})^{1/2}$
$\displaystyle<\epsilon,$

where
$\epsilon=\epsilon_{a}+\epsilon_{r}\nu_{1,1}^{-1/2}\xi_{\delta_{1}}^{\mathcal{L}}(\Delta_{1,1};x_{1,1},\nu_{1,1})^{1/2}$,
for some $\epsilon_{a}\geq 0$. In our experiments, we chose
$\epsilon_{a}=\epsilon_{r}=10^{-4}$. We could not base our criteria upon
$\xi^{\mathcal{L}}$ in (13) because we do not know the $\nu$ associated to
this measure (which is different from the $\nu_{k,j}$ generated by Algorithm 2
associated to the barrier subproblem).

Once Algorithm 3 terminates, we use a crossover technique to set $x_{i}=0$
or/and $z_{i}=0$, to respect the complementarity condition. To do so, we check
the final value of $x_{i}$ (resp. $z_{i}$), and if it is smaller than
$\sqrt{\mu}$ we set it to zero. If both $x_{i}$ and $z_{i}$ are smaller than
$\mu^{1/4}$, we set them to zero.

The subproblems in Algorithm 3 are solved with the algorithm R2 [2], and we
compare Algorithm 1 to TR [2] with R2 used as a subsolver, and R2 used by
itself. Algorithm 3 will be denoted RIPM. Finally, we introduce a variant of
RIPM named RIPMDH (_Regularized Interior Proximal Method with Diagonal Hessian
approximations_), that uses the same idea as our algorithm TRDH [17]: instead
of using LBFGS or LSR1 quasi-Newton approximations for $B_{k,j}$, we use
diagonal quasi-Newton approximations so that (25) with $\varphi$ defined in
(48) can be solved analytically for specific seperable regularizers $h$. TRDH
and RIPMDH use the Spectral Gradient update in all our results. We choose
either $h(x)=\lambda\|x\|_{0}$ or $h(x)=\lambda\|x\|_{1}$, where $\lambda>0$.
When $h(x)=\lambda\|x\|_{0}$, RIPM and RIPMDH denote Algorithm 1 instead of
Algorithm 3 because $h$ is not convex.

For simplicity, we described how to solve (1) with the constraint $x\geq 0$,
but RIPM and RIPMDH are actually able to handle box constraints $\ell\leq
x\leq u$. These more general constraints can be handled with minor
modifications using the barrier function

(82)
$\tilde{\phi}_{k}(x):=-\mu_{k}\sum_{i=1}^{n}\log(x_{i}-\ell_{i})-\mu_{k}\sum_{i=1}^{n}\log(u_{i}-x_{i})$

instead of $\phi_{k}$ [9, Section $13.8$]. When $\ell_{i}=-\infty$ (resp.
$u_{i}=+\infty$) for some $i\in\\{1,\ldots,n\\}$, we remove the term
$\log(x_{i}-\ell_{i})$ (resp. $\log(u_{i}-x_{i})$) from the first (resp.
second) sum in (82).

Our results report

  * •

the final $f(x)$;

  * •

the final $h(x)/\lambda$, where $\lambda$ is a parameter relative to our
regularization function $h$;

  * •

the final stationarity measure $\sqrt{\xi/\nu}$;

  * •

$\|x-x^{*}\|$, where $x_{*}$ is the exact solution, if it is available;

  * •

the number of smooth objective evaluations $\\#f$;

  * •

the number of gradient evaluation $\\#\nabla f$;

  * •

the number of proximal operator evaluations $\\#\textup{prox}$;

  * •

the elapsed time $t$ in seconds.

Our main goal is to reduce the number of objective and gradient evaluation, as
they are typically costly to evaluate. Since we did not fully optimize the
allocations in our algorithms, we do not pay attention to the elapsed time,
and we only report it in the tables for information.

Once a problem has been solved by all solvers, we compare their final
objective values and we save the smallest, that we denote $(f+h)^{*}$. Then,
for all solvers, we plot $(f+h)(x_{k})-(f+h)^{*}$ for every iteration $k$
where the gradient $\nabla f$ is evaluated. This allows us to represent the
evolution of the objective per gradient evaluation. For the last gradient
evaluation of RIPM and RIPMDH, we display their final objective value after
applying the crossover technique.

###  6.1 Box-constrained quadratic problem

For our first numerical experiment, we solve

(83)
$\mathop{\textup{minimize}}_{x}c^{T}x+\tfrac{1}{2}x^{T}Hx+h(x)\quad\mathop{\textup{subject
to}}\ \ell\leq x\leq u,$

which is similar to [23, Section $7.1$], where $h=\lambda\|\cdot\|_{1}$,
$H=A+A^{T}$, $A\in\mathds{R}^{n\times n}$ has nonzero components with
probability $p=10^{-4}$ following a normal law of mean $0$ and standard
deviation $1$, $c\in\mathds{R}^{n}$ has components generated using a normal
distribution of mean $0$ and standard deviation $1$, $\ell=-e-t_{\ell}$ and
$u=e+t_{u}$, with $t_{\ell}\in\mathds{R}^{n}$, $t_{u}\in\mathds{R}^{n}$ are
vectors sampled from a uniform distribution between $0$ and $1$. We chose
$n=10^{5}$, and use a LSR1 quasi-Newton approximation for TR and RIPM. For
$\lambda\geq 1.0$, the components $x_{i}$ of the solutions returned by TR,
TRDH, R2, RIPM and RIPMDH satisfy $x_{i}\in\\{\ell_{i},u_{i},0\\}$ for almost
all $i\in\\{1,\ldots,n\\}$. In this case, we observe that TR, TRDH and R2 are
more efficient than RIPM. However, as we decrease $\lambda$, we get more
components $x_{i}\notin\\{\ell_{i},u_{i},0\\}$. We show results with
$\lambda=10^{-1}$ in Figure 1 and Table 1. TRDH performs the least amount of
objective, gradient and proximal operators evaluations. RIPMDH finds the
smallest final objective value, and performs fewer objective, gradient and
proximal operator evaluations than TR-R2. RIPM terminates with a criticality
measure higher than the other solvers, but we observe that its final objective
is smaller than those of TR, TRDH and R2. For RIPM and RIPMDH, we can clearly
see plateaus that delimit the outer iterations. RIPM-R2 performs many more
proximal operator evaluations than RIPMDH, because it uses up to $200$ R2
iterations to solve (25). The number of proximal operator evaluations with
RIPM-R2 is also much higher than that of TR-R2, because the subproblems solved
with R2 in RIPM-R2 have their objective based upon (48), which is not well
conditioned when some components of $x_{k,j}$ approach $0$, whereas the
subproblems in TR-R2 are based upon (26).

Table 1:  Statistics of (83). TR and RIPM use an LSR1 Hessian approximation. The maximum number of objective evaluations is set to $800$.  solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\\#f$ | $\\#\nabla f$ | $\\#prox$ |  $t$ ($s$)
---|---|---|---|---|---|---|---
R2 |  $-2.29$e$+04$ |  $1.5$e$+04$ |  $8.5$e$-03$ | $679$ | $520$ | $679$ |  $6.8$e$-01$
TRDH |  $-2.28$e$+04$ |  $1.5$e$+04$ |  $5.9$e$-05$ | $57$ | $47$ | $113$ |  $3.6$e$-01$
TR-R2 |  $-2.28$e$+04$ |  $1.5$e$+04$ |  $9.9$e$-03$ | $801$ | $596$ | $12639$ |  $8.5$e$+00$
RIPM-R2 |  $-2.30$e$+04$ |  $1.4$e$+04$ |  $3.4$e$+00$ | $801$ | $628$ | $101019$ |  $5.0$e$+01$
RIPMDH |  $-2.32$e$+04$ |  $1.5$e$+04$ |  $8.7$e$-03$ | $313$ | $241$ | $628$ |  $2.7$e$+00$
Figure 1:  Plots of the objective of (83) per gradient evaluation with
different solvers.

###  6.2 Sparse nonnegative matrix factorization (NNMF)

The second experiment considered is the sparse nonnegative matrix
factorization (NNMF) problem from Kim and Park [15]. Let
$A\in\mathds{R}^{m\times n}$ have nonnegative entries. Each column of $A$
represents an observation, and is generated using a mixture of Gaussians where
negative entries are set to zero. We factorize $A\approx WH$ by separating $A$
into $k<\min(m,n)$ clusters, where $W\in\mathds{R}^{m\times k}$,
$H\in\mathds{R}^{k\times n}$ both have nonnegative entries and $H$ is sparse.
This problem can be written as

(84)
$\mathop{\textup{minimize}}_{W,H}\tfrac{1}{2}\|A-WH\|_{F}^{2}+h(H)\quad\mathop{\textup{subject
to}}\ W,H\geq 0,$

where $h(H)=\lambda\|\textup{vec}(H)\|_{1}$ and $\textup{vec}(H)$ stacks the
columns of $H$ to form a vector.

We set $m=100$, $n=50$, $k=5$, $\lambda=10^{-1}$, and report the statistics in
Table 2. For this particular problem, we use $\epsilon_{r}=10^{-6}$, which
allows for more accurate solves and for a better visualization of the
evolution of the objective values, shown in Figure 2. We observe that RIPM-R2
and RIPMDH are the only solvers to terminate. They outperform R2 and TR-R2 in
terms of number of objective and gradient evaluations, and their final
objective value is also smaller. R2 and TR-R2 reach the maximum number of
iterations. The objective of RIPM-R2 and RIPMDH is higher than that of R2 and
TR-R2 only in the early iterations, because the barrier function has more
effect when $\mu_{k}$ is larger. RIPM-R2 performs less objective and gradient
evaluations than RIPMDH, but much more proximal operator evaluations because
of the reasons evoked in Section 6.1.

Table 2:  Statistics of (84). TR and RIPM use an LSR1 Hessian approximation. The maximum number of objective evaluations is set to $8000$. solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\\#f$ | $\\#\nabla f$ | $\\#prox$ |  $t$ ($s$)
---|---|---|---|---|---|---|---
TRDH |  $1.25$e$+02$ |  $3.1$e$+01$ |  $1.6$e$-01$ | $8001$ | $6156$ | $16000$ |  $2.0$e$+00$
TR-R2 |  $1.25$e$+02$ |  $2.8$e$+01$ |  $1.2$e$-01$ | $8001$ | $5122$ | $150563$ |  $6.4$e$+00$
RIPM-R2 |  $1.25$e$+02$ |  $1.9$e$+01$ |  $2.5$e$-02$ | $4501$ | $3210$ | $470975$ |  $1.1$e$+01$
RIPMDH |  $1.25$e$+02$ |  $2.0$e$+01$ |  $1.4$e$-02$ | $4602$ | $3759$ | $9205$ |  $1.3$e$+00$
Figure 2:  Plots of the objective of (84) per gradient evaluation with
different solvers.

###  6.3 FitzHugh-Nagumo problem (FH)

We sample the functions $V(t;x)$ and $W(t;x)$ satisfying the FitzHugh [13] and
Nagumo et al. [20] model for neuron activation, where $x\in\mathds{R}^{5}$, as
$v(x)=(v_{1}(x),\ldots,v_{n+1}(x))$ and $w(x)=(w_{1}(x),\ldots,w_{n+1}(x))$.

(85)
$\frac{\mathrm{d}V}{\mathrm{d}t}=(V-V^{3}/3-W+x_{1})x_{2}^{-1},\quad\frac{\mathrm{d}W}{\mathrm{d}t}=x_{2}(x_{3}V-x_{4}W+x_{5}).$

The time interval $t\in[0,20]$ is discretized, with the initial conditions
$(V(0),W(0))=(2,0)$. We solve

(86)
$\mathop{\textup{minimize}}_{x}\tfrac{1}{2}\|(v(x)-\bar{v}(\bar{x}),w(x)-\bar{w}(\bar{x}))\|_{2}^{2}+h(x),\quad\mathop{\textup{subject
to}}\ x_{2}\geq 0.5,$

where $h(x)=\lambda\|x\|_{0}$ with $\lambda=10$, $n=100$, and report the
statistics in Table 3. Since $h$ is not convex, we use Algorithm 1 instead of
Algorithm 3. We do not show results with R2 because it encounters a numerical
error during the solve of a differential equation to compute the objective.
The evolution of the objective per gradient evaluation is shown in Figure 3.
To improve readability, we choose to show the number of gradient evaluations
on a logarithmic scale, and not to plot results with TRDH. All solvers
converge to the value $(0.00,0.50,0.54,0.00,0.00)$ except for TRDH that has a
higher final objective value than the other solvers. TR-R2 is the fastest, and
seems the most suited to solve smaller problems such as (86). RIPM and RIPMDH
still converge, but the latter is much slower. However, RIPMDH performs the
least amount of proximal operator evaluations.

Table 3:  Statistics of (86). TR and RIPM use an LBFGS Hessian approximation.  solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\\#f$ | $\\#\nabla f$ | $\\#prox$ |  $t$ ($s$)
---|---|---|---|---|---|---|---
TRDH |  $6.05$e$+00$ | $3$ |  $2.9$e$+01$ | $1001$ | $697$ | $2000$ |  $4.5$e$+00$
TR-R2 |  $4.40$e$+00$ | $2$ |  $4.8$e$-03$ | $53$ | $45$ | $4627$ |  $3.3$e$-01$
RIPM-R2 |  $4.40$e$+00$ | $2$ |  $1.3$e$-02$ | $261$ | $112$ | $10139$ |  $1.0$e$+00$
RIPMDH |  $4.40$e$+00$ | $2$ |  $1.5$e$-02$ | $798$ | $523$ | $1600$ |  $3.8$e$+00$
Figure 3:  Plots of the objective of (86) per gradient evaluation with
different solvers.

###  6.4 Constrained basis pursuit denoise (BPDN)

We solve the basis pursuit denoise problem (BPDN) [24, 12] with additional
bound constraints. Let $m=200$, $n=512$, $b=Ax_{\star}+\epsilon$, where
$\epsilon\sim\mathcal{N}(0,0.01)$, $A\in\mathds{R}^{m\times n}$ has
orthonormal rows, and $x_{\star}$ is a vector of zeros, except for $5$ of its
components that are set to $1$. The constrained BPDN problems is written as

(87)
$\mathop{\textup{minimize}}_{x}\tfrac{1}{2}\|Ax-b\|_{2}^{2}+h(x)\quad\mathop{\textup{subject
to}}\ x\geq 0,$

where $h(x)=\lambda\|x\|_{1}$. We use $\lambda=\|A^{T}b\|_{\infty}/10$.

The statistics are shown in Table 4. R2, TRDH and TR-R2 are much more
efficient than RIPM on this problem. This could come from the fact that there
are many active bounds in the solution. However, this was also the case for
the NNMF problem of Section 6.2, for which RIPM seems more efficient. Further
investigations should seek to understand such behaviours on different
problems. RIPM-R2-p and RIPMDH-p use the modifications $\mu_{0}=10^{-3}$ and
$\epsilon_{r,i}=1.0$, which make RIPMDH on (87) surpass TR-R2 and close to
TRDH. Figure 4 shows the evolution of the objective values. RIPM and RIPMDH
are not included to improve readability.

Table 4:  Statistics of (87). TR and RIPM use an LSR1 Hessian approximation.  solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\|x-x^{*}\|_{2}$ | $\\#f$ | $\\#\nabla f$ | $\\#prox$ |  $t$ ($s$)
---|---|---|---|---|---|---|---|---
R2 |  $3.91$e$-02$ |  $8.7$e$+00$ |  $1.8$e$-04$ |  $4.1$e$-01$ | $11$ | $11$ | $11$ |  $4.0$e$-03$
TRDH |  $3.91$e$-02$ |  $8.7$e$+00$ |  $6.8$e$-05$ |  $4.1$e$-01$ | $9$ | $9$ | $17$ |  $8.0$e$-03$
TR-R2 |  $3.91$e$-02$ |  $8.7$e$+00$ |  $1.6$e$-04$ |  $4.1$e$-01$ | $17$ | $17$ | $35$ |  $1.1$e$-02$
RIPM-R2 |  $4.12$e$-02$ |  $8.7$e$+00$ |  $7.1$e$-03$ |  $4.2$e$-01$ | $915$ | $915$ | $12143$ |  $3.3$e$+00$
RIPMDH |  $3.92$e$-02$ |  $8.7$e$+00$ |  $9.9$e$-04$ |  $4.1$e$-01$ | $361$ | $225$ | $724$ |  $2.1$e$-01$
RIPM-R2-p |  $4.22$e$-02$ |  $8.7$e$+00$ |  $1.0$e$-02$ |  $4.3$e$-01$ | $56$ | $56$ | $8331$ |  $1.5$e$-01$
RIPMDH-p |  $3.91$e$-02$ |  $8.7$e$+00$ |  $5.6$e$-04$ |  $4.1$e$-01$ | $15$ | $15$ | $32$ |  $9.0$e$-03$
Figure 4:  Plots of the objective of (87) per gradient evaluation with
different solvers.

##  7 Discussion and future work

We have presented RIPM, a trust-region interior-point method to solve
nonsmooth regularized problems with box constraints, and RIPMDH, a variant
based upon techniques of Leconte and Orban [17]. These algorithms solve a
sequence of unconstrained barrier subproblems to obtain a sequence of
approximate solutions of (1). We have shown the convergence of the inner
barrier subproblems, and we have characterized the degree of
$(\epsilon_{p},\epsilon_{d})$-KKT optimality for every outer iteration past a
certain rank. Under the assumption that the iterates remain bounded, we have
shown that RIPM converges to a first-order stationary point for (1). We
compared RIPM and RIPMDH to projected-direction methods with a separable
regularization function.

RIPM and RIPMDH perform well on the box-constrained quadratic problem of
Section 6.1 and on the NNMF problem of Section 6.2. They are not as efficient
on the FH problem of Section 6.3 and the constrained BPDN problem of Section
6.4, which may suggest that projected-direction methods may be more efficient
to solve problems with fewer variables and constraints. However, as observed
with RIPM-R2-p and RIPMDH-p, the modification of two parameters of RIPM and
RIPMDH improves their efficiency significantly on the constrained BPDN
problem. This suggests that our implementation could benefit from parameter
tuning.

Future work may include generalizing the algorithm to constraints of the form
$c_{i}(x)\leq 0$ with $i\in\\{1,\ldots,m\\}$ for some $m>0$, where the $c_{i}$
are continuously differentiable and Lipschitz-gradient continuous, as in [9,
Section $13.9$] in the smooth case, or [11] for nonsmooth problems.

Another improvement would be to scale the trust region to allow greater search
directions along the boundary of the feasible domain. This is explained more
in detail in [9, Section $13.7$] for trust-regions based upon the
$\ell_{2}$-norm. However, we could not find an alternative for trust-regions
based upon the $\ell_{\infty}$-norm that led to satisfying numerical results.

In Section 5.4, Model Assumption 5.5 does not allow the use of
$h=\|\cdot\|_{0}$ with Algorithm 3. It would be interesting to see whether it
is possible to establish convergence properties similar to those of Algorithm
1 without this assumption. One way to do this might be to replace
$\xi_{\textup{cp}}$ by $\xi_{\delta_{k}}^{\mathcal{L}}$ in (35b) when $j$ is
large enough, but we did not manage to justify that this change results in a
reasonable Step Assumption 5.1.

The extension of the convergence results to locally Lipschitz-gradient
continuous functions $f$ could also be studied, based upon the work of [10,
11, 14].

Finally, when $\\{B_{k,j}\\}_{j}$ grows unbounded, we may still be able to
prove the convergence of RIPM using the analysis of Leconte and Orban [18],
provided that the norm of the Hessian approximations do not grow too fast.

## References

  * Aravkin et al. [2023] A. Aravkin, R. Baraldi, G. Leconte, and D. Orban.  Corrigendum: A proximal quasi-Newton trust-region method for nonsmooth regularized optimization.  Les Cahiers du GERAD G-2021-12-SM, Groupe d’études et de recherche en analyse des décisions, GERAD, Montréal QC H3T 2A7, Canada, Aug. 2023.
  * Aravkin et al. [2022] A. Y. Aravkin, R. Baraldi, and D. Orban.  A proximal quasi-Newton trust-region method for nonsmooth regularized optimization.  _SIAM J. Optim._ , 32(2):900–929, 2022.
  * Attouch [1977] H. Attouch.  Convergence de fonctions convexes, des sous-différentiels et semi-groupes associés.  _CR Acad. Sci. Paris_ , 284(539-542):13, 1977\.
  * Attouch and Wets [1981] H. Attouch and R. J.-B. Wets.  Approximation and convergence in nonlinear optimization.  In O. L. Mangasarian, R. R. Meyer, and S. M. Robinson, editors, _Nonlinear Programming 4_ , pages 367–394. Academic Press, 1981.  ISBN 978-0-12-468662-5.
  * Baraldi et al. [2022] R. Baraldi, G. Leconte, and D. Orban.  RegularizedOptimization.jl: Algorithms for regularized optimization.  https://github.com/JuliaSmoothOptimizers/RegularizedOptimization.jl, February 2022.
  * Bertocchi et al. [2020] C. Bertocchi, E. Chouzenoux, M.-C. Corbineau, J.-C. Pesquet, and M. Prato.  Deep unfolding of a proximal interior point method for image restoration.  _Inverse Problems_ , 36(3):034005, feb 2020.
  * [7] E. Chouzenoux, M.-C. Corbineau, and J.-C. Pesquet.  A Proximal Interior Point Algorithm with Applications to Image Processing.  _J. Math. Imaging Vis._ , 62(6):919–940.
  * Conn et al. [2000a] A. R. Conn, N. I. M. Gould, D. Orban, and P. L. Toint.  A primal-dual trust-region algorithm for non-convex nonlinear programming.  _Math. Program._ , 87(2):215–249, 2000a.
  * Conn et al. [2000b] A. R. Conn, N. I. M. Gould, and Ph. L. Toint.  _Trust-Region Methods_.  Number 1 in MOS-SIAM Series on Optimization. SIAM, Philadelphia, USA, 2000b.
  * [10] A. De Marchi and A. Themelis.  Proximal Gradient Algorithms Under Local Lipschitz Gradient Continuity.  _J. Optim. Theory and Applics._ , 194(3):771–794.
  * De Marchi and Themelis [2024] A. De Marchi and A. Themelis.  An interior proximal gradient method for nonconvex optimization.  arXiv preprint 2208.00799, 2024.
  * Donoho [2006] D. Donoho.  Compressed sensing.  _IEEE T. Inform. Theory_ , 52(4):1289–1306, 2006\.
  * FitzHugh [1955] R. FitzHugh.  Mathematical models of threshold phenomena in the nerve membrane.  _B. Math. Biophys._ , 17(4):257–278, 1955.
  * Kanzow and Mehlitz [2022] C. Kanzow and P. Mehlitz.  Convergence Properties of Monotone and Nonmonotone Proximal Gradient Methods Revisited.  _J. Optim. Theory and Applics._ , 195(2):624–646, 2022.
  * Kim and Park [2008] J. Kim and H. Park.  Sparse nonnegative matrix factorization for clustering.  Technical Report GT-CSE-08-01, Georgia Inst. of Technology, 2008.
  * Le [2013] H. Y. Le.  Generalized subdifferentials of the rank function.  _Optimization Letters_ , 7(4):731–743, 2013.
  * Leconte and Orban [2023a] G. Leconte and D. Orban.  The indefinite proximal gradient method.  Les Cahiers du GERAD G-2023-37, Groupe d’études et de recherche en analyse des décisions, GERAD, Montréal QC H3T 2A7, Canada, Aug. 2023a.
  * Leconte and Orban [2023b] G. Leconte and D. Orban.  Complexity of trust-region methods with unbounded Hessian approximations for smooth and nonsmooth optimization.  Les Cahiers du GERAD G-2023-65, Groupe d’études et de recherche en analyse des décisions, GERAD, Montréal QC H3T 2A7, Canada, Dec. 2023b.
  * Lions and Mercier [1979] P.-L. Lions and B. Mercier.  Splitting algorithms for the sum of two nonlinear operators.  _SIAM J. Numer. Anal._ , 16(6):964–979, 1979\.
  * Nagumo et al. [1962] J. Nagumo, S. Arimoto, and S. Yoshizawa.  An active pulse transmission line simulating nerve axon.  _Proceedings of the IRE_ , 50(10):2061–2070, 1962\.
  * Poliquin [1992] R. A. Poliquin.  An extension of Attouch’s theorem and its application to second-order epi-differentiation of convexly composite functions.  _Transactions of the American Mathematical Society_ , 332:861–874, 1992.
  * Rockafellar and Wets [1998] R. Rockafellar and R. Wets.  _Variational Analysis_, volume 317.  Springer Verlag, 1998.
  * Shen et al. [2020] C. Shen, W. Xue, L.-H. Zhang, and B. Wang.  An Active-Set Proximal-Newton Algorithm for $\ell_{1}$ Regularized Optimization Problems with Box Constraints.  _Journal of Scientific Computing_ , 85(3):57, 2020\.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _J. Roy. Statist. Soc. Ser. B_ , 58(1):267–288, 1996.