# The Indefinite Proximal Gradient Method

Geoffroy Leconte111GERAD and Department of Mathematics and Industrial
Engineering, Polytechnique Montréal. E-mail<EMAIL_ADDRESS>Dominique Orban222GERAD and Department of Mathematics and Industrial
Engineering, Polytechnique Montréal. E-mail<EMAIL_ADDRESS>Research supported by an NSERC Discovery grant.

###### Abstract

We introduce a variant of the proximal gradient method in which the quadratic
term is diagonal but may be indefinite, and is safeguarded by a trust region.
Our method is a special case of the proximal quasi-Newton trust-region method
of Aravkin et al. [3]. We provide closed-form solution of the step computation
in certain cases where the nonsmooth term is separable and the trust region is
defined in the infinity norm, so that no iterative subproblem solver is
required. Our analysis expands upon that of [3] by generalizing the trust-
region approach to problems with bound constraints. We provide an efficient
open-source implementation of our method, named TRDH, in the Julia language in
which Hessians approximations are given by diagonal quasi-Newton updates. TRDH
evaluates one standard proximal operator and one indefinite proximal operator
per iteration. We also analyze and implement a variant named iTRDH that
performs a single indefinite proximal operator evaluation per iteration. We
establish that iTRDH enjoys the same asymptotic worst-case iteration
complexity as TRDH. We report numerical experience on unconstrained and bound-
constrained problems, where TRDH and iTRDH are used both as standalone and
subproblem solvers. Our results illustrate that, as standalone solvers, TRDH
and iTRDH improve upon the quadratic regularization method R2 of [3] but also
sometimes upon their quasi-Newton trust-region method, referred to here as
TR-R2, in terms of smooth objective value and gradient evaluations. On
challenging nonnegative matrix factorization, binary classification and data
fitting problems, TRDH and iTRDH used as subproblem solvers inside TR improve
upon TR-R2 for at least one choice of diagonal approximation.

##  1 Introduction

We consider the nonsmooth regularized problem

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

where $\ell\in(\mathds{R}\cup\\{-\infty\\})^{n}$,
$u\in(\mathds{R}\cup\\{+\infty\\})^{n}$ with $\ell\leq u$ componentwise,
$f:\mathds{R}^{n}\to\mathds{R}$ is continuously differentiable on an open set
containing $[\ell,u]$, and $h:\mathds{R}^{n}\to\mathds{R}\cup\\{+\infty\\}$ is
proper and lower semicontinuous (lsc). A component $\ell_{i}=-\infty$ or
$u_{i}=+\infty$ indicates that $x_{i}$ is unbounded below or above,
respectively. Both $f$ and $h$ may be nonconvex. Typically, $h$ is nonsmooth
and serves to identify a local minimizer of $f$ with desirable features, such
as sparsity.

Numerical methods for (1) are typically based on the proximal-gradient method
[21]. Aravkin et al. [3] provide an overview of recent works focusing on (1)
where both $f$ and $h$ may be nonconvex, to which we refer the reader. In
addition, they propose two methods: an adaptive quadratic regularization
approach named R2, which may be viewed as a proximal-gradient method with
adaptive step size, and a quasi-Newton trust-region method named TR in which
the subproblem consist in minimizing a quadratic approximation of $f$ about
the current iterate regularized by a model of $h$ inside a trust region.
Typically, an explicit solution to the subproblem is not known, and an
iterative procedure must be used—the authors use R2. Aravkin et al. [2]
develop a similar approach designed for applications where $f$ is a least-
squares residual. They propose a trust-region and a regularization approach.
Again, they use R2 as subproblem solver.

In both R2 and the traditional proximal-gradient method, a uniformly
positively scaled $\nabla f$ is used to compute a step, and that computation
relies on the proximal operator—see Section 2 for precise definitions. Some
authors consider positive-definite diagonal gradient scaling—see below. In the
present research, we contend that generalizing the definition of the proximal
operator by allowing a diagonal scaling of $\nabla f$, and even permitting
negative scaling factors, continues to allow us to derive analytical solutions
for several nonsmooth terms of interest in applications. Moreover, such
potentially indefinite scaling might allow the model to better capture
inherent nonconvexity in $f$ about the current iterate. We devise a trust-
region method based on those ideas and name it TRDH, which stands for _trust-
region method with diagonal Hessian approximations_. At every iteration, our
method performs the evaluation of both a classical proximal operator and a
generalized proximal operator with indefinite diagonal scaling. However, it is
possible to modify TRDH slightly to devise a variant that requires a single
generalized proximal operator evaluation per iteration while preserving the
asymptotic worst-case evaluation complexity bound. We name the variant iTRDH,
which stands for _indefinite trust-region method with diagonal Hessian_.

Diagonal gradient scaling could be referred to as diagonal quasi-Newton, and
though the literature appears to be thin on that subject, there exist a few
references. Diagonal Hessian approximations range from a multiple of the
identity, as in the traditional proximal-gradient or the spectral gradient
method, to approximations computed based on a relaxed secant equation. We go
into more details in Section 4.

Because we use the $\ell_{\infty}$-norm to define the trust region, we are
able to treat bound constraints naturally, by taking the intersection of
$[\ell,\,u]$ with the trust region. Our generalized _indefinite_ proximal
operators must take the indicator of the resulting box into account together
with $h$, or a model thereof.

A by-product of the present research is an efficient software implementation
of TRDH and iTRDH, both as standalone solvers, and as subproblem solvers for
use inside TR. TRDH and iTRDH can use several diagonal Hessian approximations:
the multiple of the identity given by the spectral gradient approximation, and
two diagonal quasi-Newton approximations based on the weak secant equation.
Our open-source implementations in Julia are available from
https://github.com/JuliaSmoothOptimizers/RegularizedOptimization.jl.

We report numerical experience conducted with TRDH, iTRDH, TR-TRDH and TR-
iTRDH on a range of unconstrained and bound-constrained problems, where TR-
TRDH and TR-iTRDH refer to TRDH and iTRDH used as subsolvers in TR. Our
general conclusion is that for at least one choice of a diagonal Hessian
approximation, TRDH and iTRDH outperform R2 in terms of evaluations of $f$ and
$\nabla f$. In addition, for at least one choice of a diagonal Hessian
approximation, TR-TRDH and TR-iTRDH outperform TR-R2 on the same metric. Our
detailed results are in Section 6.

### Related research

Most of the literature focuses on positive-definite diagonal scaling of the
proximal operator; Becker and Fadili [9] and Becker et al. [10] consider
positive-definite quasi-Newton approximations of the form $H=D+VV^{T}$ where
$D$ is positive definite and diagonal, and a specialized procedure to solve
the proximal quasi-Newton subproblems. Under their assumptions, $H$ remains
uniformly bounded. Their analysis is restricted to $f$ and $h$ convex and does
not provide complexity bounds.

Duchi et al. [17] present the ADAGRAD algorithm, which is commonly used in
online learning. ADAGRAD is a stochastic algorithm but a deterministic
implementation of it would minimize a smooth objective such as $f$ using a
variant of the projected gradient algorithm, that can be seen as a special
case of the proximal gradient method, using a positive diagonal gradient
scaling. At iteration $k$, the scaling is set to
$\mathop{\text{diag}}(G_{k})^{-1/2}$ where $G_{k}=\sum_{j=0}^{k}\nabla
f(x_{j})\nabla f(x_{j})^{T}$.

Scheinberg and Tang [27] also focus on positive-definite and uniformly-bounded
Hessian approximations for $f$ and $h$ convex. Their numerical results employ
a limited-memory BFGS approximation, although the latter is liable to grow
unbounded [14, §$8.4$]. The proximal quasi-Newton subproblem is solved
inexactly with a coordinate descent algorithm.

Aravkin et al. [3] propose a proximal quasi-Newton trust-region method for
nonconvex $f$ and $h$ under weak assumptions, accompanied by a complexity
analysis. They employ limited-memory BFGS or SR1 approximations and solve
subproblems using the proximal-gradient method or R2. Though their analysis
assumes the Hessian approximations are uniformly bounded, there are known
procedures to estimate bounds on the eigenvalues of quasi-Newton
approximations after an update, and those bounds can be used to skip or modify
the update to keep it bounded [22].

Similarly, Baraldi and Kouri [4] present a proximal trust-region method with
inexact objective and gradient evaluations, under the additional assumption
that $h$ must be convex. They show numerical results with exact Hessian when
it is available, and employ iterative solvers to approximate it otherwise.

### Notation

$\mathds{B}$ is the unit ball centered at the origin and defined by a norm
dictated by the context, and $\Delta\mathds{B}$ is the ball of radius
$\Delta>0$ centered at the origin. For fixed $x\in\mathds{R}^{n}$, the ball of
radius $\Delta$ centered at $x$ is $x+\Delta\mathds{B}$. When it is necessary
to indicate that $\mathds{B}$ is defined by the $\ell_{p}$-norm, for $1\leq
p\leq+\infty$, we write $\mathds{B}_{p}$. For $A\subseteq\mathds{R}^{n}$, the
indicator of $A$ is $\chi(\cdot\mid
A):\mathds{R}^{n}\to\mathds{R}\cup\\{+\infty\\}$ defined as $\chi(x\mid A)=0$
if $x\in A$ and $+\infty$ otherwise. If $A\neq\varnothing$, $\chi(\cdot\mid
A)$ is proper. If $A$ is closed, $\chi(\cdot\mid A)$ is lsc. For
$D\in\mathds{R}^{m\times n}$ with elements $\delta_{i,j}$,
$|D|\in\mathds{R}^{m\times n}$ has elements $|\delta_{i,j}|$. For a finite set
$A\subset\mathds{N}$, we denote $|A|$ its cardinality. If $f_{1}$ and $f_{2}$
are two positive functions of $\epsilon>0$, we say that
$f_{1}(\epsilon)=O(f_{2}(\epsilon))$ if there exists a constant $C>0$ such
that $f_{1}(\epsilon)\leq Cf_{2}(\epsilon)$ for all $\epsilon>0$ sufficiently
small.

##  2 Background

The following are standard variational analysis concepts—see, e.g., [26]. 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 $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 $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$.
$\widehat{\partial}\phi(\bar{x})\subset\partial\phi(\bar{x})$ always holds.

The horizon subdifferential of $\phi$ at $\bar{x}$ is the closed, but not
necessarily convex, cone $\partial^{\infty}\phi(\bar{x})$ of
$v\in\mathds{R}^{n}$ for which there exist $\\{x_{k}\\}\to\bar{x}$,
$\\{v_{k}\\}$ and $\\{\lambda_{k}\\}\downarrow 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 $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 called
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 called 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 $\phi:\mathds{R}^{n}\to\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}$
is proper lsc, and $C\subseteq\mathds{R}^{n}$ is closed, we say that the
_constraint qualification_ is satisfied at $\bar{x}\in C$ for the constrained
problem

$\mathop{\textup{minimize}}_{x\in\mathds{R}^{n}}\ \phi(x)\quad\textup{subject
to}\ x\in C$

if

$\partial^{\infty}\phi(\bar{x})\cap N_{C}(\bar{x})=\\{0\\}.$ (2)

As an example where (2) fails to hold, let $n=1$, $f(x)=0$ for all
$x\in\mathds{R}$,

$h(x)=\begin{cases}x&\text{if }x\geq 0,\\\ +\infty&\text{if }x<0,\end{cases}$

and $C=[0,\,1]$. Clearly, $\bar{x}=0$ is the only stationary point and
$N_{C}(\bar{x})=\\{v\mid v\leq 0\\}$. Let $\mathop{\textup{epi}}h$ denote the
epigraph of $h$, i.e., the set $\\{(x,t)\mid x\in\mathds{R}^{n},\ t\geq
h(x)\\}$. Because $h$ is locally lsc about $\bar{x}$, [26, Theorem $8.9$]
yields

$\partial^{\infty}h(\bar{x})=\\{v\mid(v,0)\in
N_{\mathop{\textup{epi}}h}(\bar{x},h(\bar{x}))\\}=\\{v\mid(v,0)\in
N_{\mathop{\textup{epi}}h}(0,0)\\}=N_{C}(\bar{x}),$

and the constraint qualification does not hold. The reason is that $h$ and
$\chi(\cdot\mid C)$ both have a jump discontinuity at $\bar{x}$ and take the
value $+\infty$ for $x<\bar{x}$, which is akin to repeating a constraint.

We say that $\bar{x}$ is first-order stationary for (1) if
$0\in\partial(f+h+\chi(\cdot\mid[\ell,u]))(\bar{x})=\nabla
f(\bar{x})+\partial(h+\chi(\cdot\mid[\ell,u]))(\bar{x})$. If $\bar{x}$ is a
local solution of (1), it is first-order stationary [26, Theorem $10.1$].
Under (2) applied to (1), which reads $(\nabla
f(\bar{x})+\partial^{\infty}h(\bar{x}))\cap N_{[\ell,u]}(\bar{x})=\\{0\\}$,
the necessary optimality condition can be written equivalently as $0\in\nabla
f(\bar{x})+\partial h(\bar{x})+N_{[\ell,u]}(\bar{x})$ [26, Theorem $8.15$].
When (1) is unconstrained or $\bar{x}\in\mathop{\textup{int}}[\ell,u]$, the
normal cone is $\\{0\\}$, the constraint qualification is satisfied, and the
necessary condition reduces to $0\in\nabla f(\bar{x})+\partial h(\bar{x})$.

The proximal operator associated with $h$ is

$\mathop{\textup{prox}}_{\nu h}(q):=\mathop{\textup{argmin}}_{x}\
\tfrac{1}{2}\nu^{-1}\|x-q\|_{2}^{2}+h(x),$ (3)

where $\nu>0$ is a preset steplength. If $h$ is prox-bounded and $\nu>0$ is
sufficiently small, $\mathop{\textup{prox}}_{\nu h}(q)$ is a nonempty and
closed set. It may contain multiple elements.

The proximal gradient method [21] is a generalization of the gradient method
that takes the nonsmooth term into account. It generates iterates
$\\{x_{k}\\}$ according to

$x_{k+1}\in\mathop{\textup{prox}}_{\nu_{k}h}(x_{k}-\nu_{k}\nabla f(x_{k})).$
(4)

Equivalently, and more instructively,

$\displaystyle s_{k}$ $\displaystyle\in\mathop{\textup{argmin}}_{s}\
\tfrac{1}{2}\nu_{k}^{-1}\|s+\nu_{k}\nabla f(x_{k})\|_{2}^{2}+h(x_{k}+s),$ (5a)
$\displaystyle=\mathop{\textup{argmin}}_{s}\ \nabla
f(x_{k})^{T}s+\tfrac{1}{2}\nu_{k}^{-1}\|s\|_{2}^{2}+h(x_{k}+s),$ (5b)
$\displaystyle x_{k+1}$ $\displaystyle=x_{k}+s_{k},$ (5c)

where it becomes clear that the computation of $s_{k}$ results from
formulating a simple quadratic model

$f(x_{k})+\nabla f(x_{k})^{T}s+\tfrac{1}{2}\nu_{k}^{-1}\|s\|_{2}^{2}\approx
f(x_{k}+s)$ (6)

of the smooth term $f$, to which we add the nonsmooth term $h$.

The effectiveness of (5) comes from the observation that it is possible to
determine a closed form solution $s_{k}$, or at least a specialized procedure
to identify such a global minimizer, for numerous choices of $h$ that are
relevant in applications [8], without resort to a general-purpose optimization
solver.

Aravkin et al. [3] devise a trust-region quasi-Newton method for (1) without
explicit bound constraints in which subproblems are solved inexactly by a
method closely related to (5). However, the subproblem has the additional
trust-region constraint $s\in\Delta_{k}\mathds{B}$:

$s_{k}\in\mathop{\textup{argmin}}_{s}\ \nabla
f(x_{k})^{T}s+\tfrac{1}{2}\nu_{k}^{-1}\|s\|_{2}^{2}+h(x_{k}+s)+\chi(s\mid\Delta_{k}\mathds{B}),$
(7)

where $\Delta_{k}>0$, which may be construed as changing $h(x_{k}+s)$ to
$h(x_{k}+s)+\chi(s\mid\Delta_{k}\mathds{B})$. Because both $s\mapsto
h(x_{k}+s)$ and $s\mapsto\chi(s\mid\Delta\mathds{B})$ are proper and lsc, so
is their sum. In addition, because $\Delta_{k}\mathds{B}$ is bounded, $h$ no
longer needs to be prox-bounded for a solution to exist. Again, a closed form
solution may be found for a variety of cases of interest in applications [3,
2].

##  3 The indefinite proximal operator

###  3.1 Definition and Examples

We restrict ourselves to the case where $h$ is separable, i.e.,
$h(x)=h_{1}(x_{1})+\cdots+h_{n}(x_{n})$, in which case the search for an
element $s$ of (6) decouples componentwise:

$s_{i}\in\mathop{\textup{argmin}}_{s_{i}}\
g_{i}s_{i}+\tfrac{1}{2}\nu^{-1}s_{i}^{2}+h_{i}(x_{i}+s_{i}),$

where $g_{i}$ is the $i$-th component of $\nabla f(x_{k})$. If $\mathds{B}$ is
defined in the $\ell_{\infty}$-norm, (7) also decouples. Without loss of
generality, and accounting for the bound constraints of (1), we write it as

$s_{i}\in\mathop{\textup{argmin}}_{s_{i}}\
g_{i}s_{i}+\tfrac{1}{2}\nu^{-1}s_{i}^{2}+h_{i}(x_{i}+s_{i})+\chi(s_{i}\mid[\tilde{\ell}_{i},\,\tilde{u}_{i}]),$
(8)

where $\tilde{\ell}_{i}:=\max(\ell_{i}-x_{i},-\Delta)$ and
$\tilde{u}_{i}:=\min(u_{i}-x_{i},\Delta)$ are the $i$th components of
$\tilde{\ell}$ and $\tilde{u}\in\mathds{R}^{n}$, which describe the
intersection of the box $[\ell,u]$ with the trust region.

Replacing the quadratic coefficient in (8) by a positive scalar $\delta_{i}>0$
does not complicate the problem as it simply amounts to changing the value of
$\nu$. In the following, we contend that replacing the quadratic coefficient
by $\delta_{i}\leq 0$ continues to allow us to derive a closed form solution
in certain cases. Observe that modifying the quadratic term as described
amounts to changing (6) to

$f(x_{k})+\nabla f(x_{k})^{T}s+\tfrac{1}{2}s^{T}D_{k}s,$

where $D_{k}=\mathop{\text{diag}}(d_{k})$, and $d_{k}\in\mathds{R}^{n}$ has
components $\delta_{k.i}$. Crucially, we wish to allow $D_{k}$ to be
indefinite to capture any nonconvexity of $f$ about $x_{k}$ to a certain
extent. While such change may make the quadratic model nonconvex, the
indicator of $[\tilde{\ell}_{i},\,\tilde{u}_{i}]$ ensures that a finite
solution exists.

######  Definition 1.

Let $g\in\mathds{R}^{n}$, $D\in\mathds{R}^{n\times n}$ be diagonal and
$C\subset\mathds{R}^{n}$ be nonempty. The indefinite proximal operator of $h$
with respect to $g$, $D$ and $C$ is

$\mathop{\textup{iprox}}_{g,D,h,C}(q):=\mathop{\textup{argmin}}_{x}\
g^{T}x+\tfrac{1}{2}x^{T}Dx+h(x)+\chi(x\mid C).$ (9)

The usual proximal operator of $\nu h$ at $q\in\mathds{R}^{n}$ is a special
case of Definition 1 in which $C=\mathds{R}^{n}$, $D=\nu^{-1}I$, and
$g=-\nu^{-1}q$.

######  Example 1.

Consider $h(x)=\lambda\|x\|_{0}$, where $\lambda>0$, and $C:=[\ell,\,u]$,
where $\ell$, $u\in\mathds{R}^{n}$ with $\ell<u$ componentwise. The
computation of $x$ in (9) decouples componentwise as

$x_{i}\in\mathop{\textup{argmin}}_{x_{i}}g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}+\lambda|x_{i}|_{0}+\chi(x_{i}\mid[\ell_{i},\,u_{i}]).$

Note that

$g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}+\lambda|x_{i}|_{0}=\begin{cases}0&\text{if
}x_{i}=0\\\ \lambda+g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}&\text{if
}x_{i}\neq 0.\end{cases}$

If $\delta_{i}=0$, the solution set of the above problem is

  * •

$[\ell_{i},u_{i}]$ if $g_{i}=0$ and either $\ell_{i}>0$ or $u_{i}<0$;

  * •

$\mathop{\textup{argmin}}\\{\lambda+g_{i}x_{i}\mid x_{i}\in M_{i}\\}$ where
$M_{i}$ is the set $\\{\ell_{i},u_{i}\\}$ together with $0$ if
$0\in(\ell_{i},u_{i})$ in all other cases.

If $\delta_{i}\neq 0$, the solution set is
$\mathop{\textup{argmin}}\\{\lambda+g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}\mid
x_{i}\in M_{i}\\}$, where $M_{i}$ is the finite set comprising $\ell_{i}$,
$u_{i}$ together with

  * •

$-g_{i}/\delta_{i}$ if $\delta_{i}>0$;

  * •

$0$ if $0\in(\ell_{i},u_{i})$.

See Figure 1 for an illustration of a few cases.

$0$$-g_{i}/\delta_{i}$

$\ell_{i}$0$u_{i}$

$0$$-g_{i}/\delta_{i}$

Figure 1:  Illustration of Example 1 with $\delta_{i}>0$ (left),
$\delta_{i}=0$ (center), and $\delta_{i}<0$ (right). The blue curve is
$g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}$ while the red curve is
$g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}+\lambda|x_{i}|_{0}$.

######  Example 2.

Consider $h(x)=\lambda\|x\|_{1}$ where $\lambda>0$, and $C$ as in Example 1.
The computation of $x$ in (9) decouples componentwise as

$x_{i}\in\mathop{\textup{argmin}}_{x_{i}}g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}+\lambda|x_{i}|+\chi(x_{i}\mid[\ell_{i},\,u_{i}]).$

Note that

$g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}+\lambda|x_{i}|=\begin{cases}(g_{i}+\lambda)x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}&\text{if
}x_{i}\geq 0\\\ (g_{i}-\lambda)x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}&\text{if
}x_{i}<0.\end{cases}$

If $\delta_{i}=0$, the solution of the above problem is

  * •

$[\ell_{i},\,\max(0,\ell_{i})]$ if $g_{i}=\lambda$;

  * •

$[\min(u_{i},0),\,u_{i}]$ if $g_{i}=-\lambda$;

  * •

$\mathop{\textup{proj}}(0\mid[\ell_{i},u_{i}])$ if $0\leq|g_{i}|<\lambda$;

  * •

$\\{\ell_{i}\\}$ if $g_{i}>\lambda$;

  * •

$\\{u_{i}\\}$ if $g_{i}<-\lambda$.

Consider now the case where $\delta_{i}>0$. The branch of quadratic for
$x_{i}\geq 0$ lies above the horizontal axis if $g_{i}\geq-\lambda$, and that
for $x_{i}<0$ does the same if $g_{i}\leq\lambda$. Thus, if
$|g_{i}|\leq\lambda$, the unique unconstrained minimizer is $0$ and, by
convexity, the solution of the constrained problem is
$\mathop{\textup{proj}}(0\mid[\ell_{i},u_{i}])$.

If $g_{i}>\lambda$, for the unconstrained minimizer is
$\bar{x}_{i}:=-(g_{i}-\lambda)/\delta_{i}<0$. By convexity, the solution of
the constrained problem is
$\mathop{\textup{proj}}(\bar{x}_{i}\mid[\ell_{i},u_{i}])$. The situation is
similar when $g_{i}<-\lambda$.

Consider finally the case where $\delta_{i}<0$. In this case, the unique
constrained minimizer is
$\mathop{\textup{argmin}}\\{g_{i}x_{i}+\tfrac{1}{2}\delta_{i}x_{i}^{2}+\lambda|x_{i}|\mid
x_{i}\in M_{i}\\}$ where $M_{i}$ is the set containing $\ell_{i}$ and $u_{i}$
together with $0$ if the latter lies inside $[\ell_{i},u_{i}]$.

Figure 2 illustrates a few representative cases.

0

0

0

0$\bar{x}_{i}$

0

$\bar{x}_{i}$0

0

0

0

Figure 2:  Illustration of Example 2. Top row: $\delta_{i}=0$. From left to
right: $0\leq g_{i}<\lambda$, $g_{i}=\lambda$, and $g_{i}>\lambda$. Middle
row: $\delta_{i}>0$. From left to right: $g_{i}<-\lambda$,
$|g_{i}|\leq\lambda$, and $g_{i}>\lambda$. Bottom row: $\delta_{i}<0$. From
left to right: $g_{i}\leq-\lambda$, $|g_{i}|<\lambda$, and $g_{i}>\lambda$.

###  3.2 Properties

At $x\in\mathds{R}^{n}$ where $h$ is finite, consider models

$\displaystyle\varphi(s;x)$ $\displaystyle\phantom{:}\approx f(x+s)$ (10a)
$\displaystyle\psi(s;x)$ $\displaystyle\phantom{:}\approx h(x+s)$ (10b)
$\displaystyle m(s;x)$ $\displaystyle:=\varphi(s;x)+\psi(s;x).$ (10c)

We make the following assumption on the models (10).

######  Model Assumption 3.1.

For any $x\in\mathds{R}^{n}$, $\varphi(\cdot;x)\in\mathcal{C}^{1}$, and
satisfies $\varphi(0;x)=f(x)$ and $\nabla\varphi(0;x)=\nabla f(x)$. For any
$x\in\mathds{R}^{n}$ where $h$ is finite, $\psi(\cdot;x)$ is proper lsc, and
satisfies $\psi(0;x)=h(x)$ and $\partial\psi(0;x)=\partial h(x)$.

######  Proposition 1.

Let $C\subset\mathds{R}^{n}$ be nonempty and compact, and let Model Assumption
3.1 be satisfied. Let (1) satisfy the constraint qualification (2) at $x\in
C$. Assume $0\in\mathop{\textup{argmin}}_{s}m(s;x)+\chi(x+s\mid C)$, and let
the latter subproblem satisfy the constraint qualification (2) at $s=0$. Then
$x$ is first-order stationary for (1).

###### Proof.

Let $\chi_{x,C}(s):=\chi(x+s\mid C)$. Because $C$ is compact, Model Assumption
3.1 ensures that $s\mapsto m(s;x)+\chi(x+s\mid C)$ is lsc, bounded below and
that a minimizer exists. By assumption, $0\in\nabla
f(x)+\partial(\psi+\chi_{x,C})(0)$. By the constraint qualification for the
subproblem and Model Assumption 3.1, $0\in\nabla f(x)+\partial h(x)+N_{C}(x)$,
which are the first-order optimality conditions for (1). ∎

In the assumptions of Proposition 1, satisfaction of the constraint
qualification in the subproblem is ensured if
$\partial^{\infty}\psi(0;x)=\partial^{\infty}h(x)$ provided (2) is satisfied
at $x$ for (1).

##  4 Diagonal quasi-Newton methods

The simplest possible diagonal approximation to $\nabla^{2}f(x_{k})$ is
$B_{k}:=\sigma_{k}I$ for some scalar $\sigma_{k}$, where $I$ denotes the
identity. In the R2 algorithm [3, Algorithm $6.1$], $\sigma_{k}$ is chosen
adaptively at each iteration based on progress. A more sophisticated strategy
consists in choosing $\sigma_{k}$ so as to best approximate
$\nabla^{2}f(x_{k})$, in a sense to be defined, based on recently observed
local information. That is the logic behind the spectral projected gradient
method [11], which is easily generalized as a spectral _proximal_ gradient
method, and where

$\sigma_{k}:=s_{k-1}^{T}y_{k-1}/s_{k-1}^{T}s_{k-1},\quad
s_{k-1}:=x_{k}-x_{k-1},\quad y_{k-1}:=\nabla f(x_{k})-\nabla f(x_{k-1}).$

This choice of $\sigma_{k}$ arises from the secant equation
$\sigma_{k}s_{k-1}=y_{k-1}$, which, in general, has no solution, but can be
solved in the least-squares sense by minimizing $\|\sigma
s_{k-1}-y_{k-1}\|_{2}$ in terms of $\sigma$.

Gilbert and Lemaréchal [19] experiment, among others, with diagonal updates of
quasi-Newton approximations. However, there may be no diagonal solution
$B_{k}$ to the secant equation $B_{k}s_{k-1}=y_{k-1}$, and one must resort to
a different approach to update a diagonal approximation. Dennis and Wolkowicz
[15] introduce the _weak secant_ equation
$s_{k-1}^{T}B_{k}s_{k-1}=s_{k-1}^{T}y_{k-1}$ and accompanying update formulae.
Nazareth [24] observes that some of those proposed updates can be constrained
to update only the diagonal of $B_{k}$ but that the resulting $B_{k}$ no
longer necessarily satisfies the weak secant equation.

Zhu et al. [30] supplement the weak secant equation by additionally requiring
that $B_{k}$ be diagonal, and call the resulting set of conditions the _quasi-
Cauchy_ conditions. In particular, the spectral gradient approximation is the
only solution of the quasi-Cauchy conditions that is a multiple of $I$. They
derive updates from variational principles resembling those used in classic
quasi-Newton updates. The first update sets $B_{k}$ to the unique solution of

$\mathop{\textup{minimize}}_{B}\ \|B-B_{k-1}\|_{F}\quad\textup{subject to}\
s_{k-1}^{T}Bs_{k-1}=s_{k-1}^{T}y_{k-1},$ (11)

which is the counterpart of the Powell-symmetric-Broyden (PSB) variational
problem. The solution of (11) is shown to be

$B_{k}=B_{k-1}+\frac{s_{k-1}^{T}(y_{k-1}-B_{k-1}s_{k-1})}{\mathop{\textup{trace}}(S_{k-1}^{4})}\,S_{k-1}^{2},\quad
S_{k-1}:=\mathop{\text{diag}}(s_{k-1}).$ (12)

Like the PSB update, (12) does not possess the hereditary positive-
definiteness property. However, capturing negative curvature is a desirable
feature in trust-region methods.

Motivated by linesearch methods, Zhu et al. [30] devise a second strategy in
which they update $B_{k}^{1/2}$, and that results in hereditary positive
definiteness. The main idea is analogous to (11) with the objective replaced
with $\|B^{1/2}-B_{k-1}^{1/2}\|_{F}$. The update depends on the root of a
scalar nonlinear equation, and is therefore more costly to perform than (12).

Andrei [1] suggests an alternative variational problem in which $B_{k}$ is
determined as the solution of

$\mathop{\textup{minimize}}_{B}\
\tfrac{1}{2}\|B-B_{k-1}\|_{F}^{2}+\mathop{\textup{trace}}(B)\quad\textup{subject
to}\ s_{k-1}^{T}Bs_{k-1}=s_{k-1}^{T}y_{k-1},$ (13)

where minimizing the trace of $B$ tends to cluster its eigenvalues. The
solution of (13) is shown to be

$B_{k}=B_{k-1}+\frac{s_{k-1}^{T}(y_{k-1}+s_{k-1}-B_{k-1}s_{k-1})}{\mathop{\textup{trace}}(S_{k-1}^{4})}\,S_{k-1}^{2}-I.$
(14)

Again, (14) does not preserve positive definiteness. Also motivated by
linesearch methods, Andrei [1] employs a thresholding strategy that ensures
computation of a descent direction.

Though Nazareth [24] appears to have been motivated by derivative-free
methods,333In the sense that the quasi-Cauchy condition only requires the
gradient via $s^{T}y$. diagonal quasi-Newton updates are good candidate
approximations for use in first-order methods for regularized optimization,
and we are not aware of any prior work using them.

##  5 The indefinite proximal gradient method

###  5.1 Description of the algorithm and convergence properties

For $\nu>0$, define

$\displaystyle\varphi_{\textup{cp}}(s;x)$ $\displaystyle:=f(x)+\nabla
f(x)^{T}s,$ (15a) $\displaystyle m(s;x,\nu)$
$\displaystyle:=\varphi_{\textup{cp}}(s;x)+\tfrac{1}{2}\nu^{-1}\|s\|_{2}^{2}+\psi(s;x).$
(15b)

Guided by Proposition 1, we compute a first step denoted $s_{k,1}$ such that

$s_{k,1}\in\mathop{\textup{argmin}}_{s}\
m(s;x_{k},\nu_{k})+\chi(x_{k}+s\mid[\ell,\,u]\cap(x_{k}+\Delta_{k}\mathds{B})),$
(16)

for an appropriate value of $\nu_{k}>0$.

Let

$m(s;x,D):=\varphi_{\textup{cp}}(s;x)+\tfrac{1}{2}s^{T}Ds+\psi(s;x),$ (17)

where $D$ is a diagonal matrix. The indefinite proximal gradient iteration
with trust region for (1) is defined by changing (16) to

$s_{k}\in\mathop{\textup{argmin}}_{s}\
m(s;x_{k},D_{k})+\chi(x_{k}+s\mid[\ell,\,u]\cap(x_{k}+\Delta_{k}\mathds{B})),$
(18)

where $D_{k}$ is diagonal, and by updating $x_{k}$ and $\Delta_{k}$ as is
customary in trust-region methods [14]. Among other possible choices, we focus
on the case where $D_{k}$ results from a diagonal quasi-Newton update.

We summarize the entire procedure as Algorithm 5.1, which is a special case of
[3, Algorithm $3.1$].

Algorithm 5.1 Nonsmooth trust-region algorithm with diagonal Hessian.

1:Choose constants  $0<\eta_{1}\leq\eta_{2}<1,\quad
0<1/\gamma_{3}\leq\gamma_{1}\leq\gamma_{2}<1<\gamma_{3}\leq\gamma_{4},\quad\text{and}\quad\alpha>0,\
\beta\geq 1.$

2:Choose $x_{0}\in\mathds{R}^{n}$ where $h$ is finite, $\Delta_{0}>0$, compute
$f(x_{0})+h(x_{0})$.

3:Choose $d_{0}\in\mathds{R}^{n}$ and set $D_{0}=\mathop{\text{diag}}(d_{0})$.

4:for $k=0,1,\dots$ do

5: Set $\nu_{k}:=1/(\|d_{k}\|_{\infty}+\alpha^{-1}\Delta_{k}^{-1})$.

6: Define $\varphi(s;x_{k})$ and $\psi(s;x_{k})$ according to Model Assumption
3.1.

7: Define $m(s;x_{k},\nu_{k})$ as in (15) and compute $s_{k,1}$ as in (16).

8: Define $m(s;x_{k},D_{k})$ as in (17) and compute a solution $s_{k}$ of (18)
with $\Delta_{k}$ replaced by $\min(\Delta_{k},\,\beta\|s_{k,1}\|)$.

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

10: If $\rho_{k}\geq\eta_{1}$, set $x_{k+1}=x_{k}+s_{k}$. Otherwise, set
$x_{k+1}=x_{k}$.

11: Choose $d_{k+1}\in\mathds{R}^{n}$ and set
$D_{k+1}=\mathop{\text{diag}}(d_{k+1})$.

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

In [3, Algorithm $3.1$], convergence hinges crucially on the fact that the
decrease in $m(\cdot;x_{k},D_{k})$ achieved by $s_{k}$ is at least a fraction
of the decrease in $m(\cdot;x_{k},\nu_{k})$ achieved by $s_{k,1}$, the first
proximal gradient step with a well-chosen step size $\nu_{k}>0$. At iteration
$k$ of Algorithm 5.1, $\nu_{k}<1/\|d_{k}\|_{\infty}$. Note that the choice
$d_{k}=0$ is allowed.

In Algorithm 5.1, the computation of $s_{k,1}$ serves two purposes. The first
is as stopping condition by Proposition 1, and the second is to set the trust-
region radius in the subproblem for $s_{k}$. This is at variance with [3,
Algorithm $3.1$], where, in addition, $s_{k}$ is computed by continuing the
proximal gradient iterations from $s_{k,1}$.

In a similar notation to that of [3], let

$\displaystyle\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})$
$\displaystyle:=f(x_{k})+h(x_{k})-(\varphi_{\textup{cp}}(s_{k,1};x_{k})+\psi(s_{k,1};x_{k})),$
(19a) $\displaystyle\xi(\Delta_{k};x_{k},D_{k})$
$\displaystyle:=f(x_{k})+h(x_{k})-m(s_{k};x_{k},D_{k})$ (19b)

denote the optimal model decrease for (15) and (17). The definition of (16)
guarantees that

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

which implies

$\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})\geq\tfrac{1}{2}\nu_{k}^{-1}\|s_{k,1}\|^{2}=\tfrac{1}{2}(\|d_{k}\|_{\infty}+\alpha^{-1}\Delta_{k}^{-1})\|s_{k,1}\|^{2}\geq\tfrac{1}{2}\alpha^{-1}\Delta_{k}^{-1}\|s_{k,1}\|^{2}.$
(21)

If $\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})=0$, we obtain $s_{k,1}=0$.
Proposition 1 then implies that $x_{k}$ is first-order stationary provided (1)
satisfies the constraint qualification at $x_{k}$ and (15) satisfies it at
$s=0$.

To establish convergence properties of Algorithm 5.1, we require Aravkin et
al.’s Step Assumption $3.1$, recalled in the following assumption for
convenience.

######  Step Assumption 5.1.

There exists $\kappa_{\textup{m}}>0$ and $\kappa_{\textup{mdc}}\in(0,1)$ such
that for all $k$,

$\displaystyle|f(x_{k}+s_{k})+h(x_{k}+s_{k})-(\varphi(s_{k};x_{k},D_{k})+\psi(s_{k};x_{k}))|$
$\displaystyle\leq\kappa_{\textup{m}}\|s_{k}\|_{2}^{2},$ (22a)
$\displaystyle\xi(\Delta_{k};x_{k},D_{k})$
$\displaystyle\geq\kappa_{\textup{mdc}}\xi_{\textup{cp}}(\Delta_{k};x_{k};\nu_{k}).$
(22b)

Among other situations, (22a) is satisfied if $\\{D_{k}\\}$ is bounded,
$\nabla f$ is Lipschitz continuous, and we select $\psi(s;x_{k})=h(x_{k}+s)$.
The following proposition gives a sufficient condition for which (22b) is
satisfied.

######  Proposition 2.

If

$\|d_{k}\|_{\infty}\leq(\kappa_{\textup{mdc}}^{-1}-1)\alpha^{-1}\Delta_{k}^{-1},$
(23)

then (22b) is satisfied.

###### Proof.

By definition of $s_{k}$,

$\displaystyle m(s_{k};x_{k},D_{k})$ $\displaystyle\leq
m(s_{k,1};x_{k},D_{k})$
$\displaystyle=\varphi(s_{k,1};x_{k})+\tfrac{1}{2}s_{k,1}^{T}D_{k}s_{k,1}+\psi(s_{k,1};x_{k})$
$\displaystyle\leq\varphi(s_{k,1};x_{k})+\tfrac{1}{2}\|d_{k}\|_{\infty}\|s_{k,1}\|^{2}+\psi(s_{k,1};x_{k})$
$\displaystyle=\varphi(s_{k,1};x_{k})+\tfrac{1}{2}(\nu_{k}^{-1}-\alpha^{-1}\Delta_{k}^{-1})\|s_{k,1}\|^{2}+\psi(s_{k,1};x_{k}),$

which leads to

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

For (22b) to be satisfied, it is sufficient to show that

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

that we rewrite as

$(1-\kappa_{\textup{mdc}})\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})\geq\tfrac{1}{2}(\nu_{k}^{-1}-\alpha^{-1}\Delta_{k}^{-1})\|s_{k,1}\|^{2}.$
(24)

Let
$\|d_{k}\|_{\infty}\leq(\kappa_{\textup{mdc}}^{-1}-1)\alpha^{-1}\Delta_{k}^{-1}$,
which is equivalent to
$\nu_{k}=\dfrac{1}{\|d_{k}\|_{\infty}+\alpha^{-1}\Delta_{k}^{-1}}\geq\kappa_{\textup{mdc}}\alpha\Delta_{k}$.
Since $\kappa_{\textup{mdc}}\in(0,1)$, we have

$(1-\kappa_{\textup{mdc}})\tfrac{1}{2}\nu_{k}^{-1}\geq\tfrac{1}{2}(\nu_{k}^{-1}-\alpha^{-1}\Delta_{k}^{-1}).$

By multiplying by $\|s_{k,1}\|^{2}$ and recalling that
$\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})\geq\tfrac{1}{2}\nu_{k}^{-1}\|s_{k,1}\|^{2}$,
we get

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

so that (24) is satisfied. ∎

We can always choose $\kappa_{\textup{mdc}}$ small enough to ensure that
$\|d_{k}\|_{\infty}$ does not have to be too close to zero if $\Delta_{k}$
does not grow unbounded. Whenever Algorithm 5.1 generates infinitely many very
successful iterations, we could define a $\Delta_{\max}$ from which we stop
increasing $\Delta_{k}$. We point out that Proposition 2 only gives a
sufficient condition for (22b) to be satisfied, but it is not necessary.
Therefore, for the rest of this paper, we base our analysis on (22b) rather
that on (23).

Under Step Assumption 5.1, we can apply directly Aravkin et al.’s convergence
properties.

######  Proposition 3 (3, Theorem $3.4$).

Let Step Assumption 5.1 be satisfied and

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

If $x_{k}$ is not first-order stationary and
$\Delta_{k}\leq\Delta_{\textup{succ}}$, iteration $k$ is very successful and
$\Delta_{k+1}>\Delta_{k}$.

Without further assumptions, the following result also holds.

######  Proposition 4 (3, Theorem $3.5$).

Let Step Assumption 5.1 be satisfied and assume Algorithm 5.1 generates a
finite number of successful iterations. Then $x_{k}=x_{\star}$ for all
sufficiently large $k$. If (1) satisfies the constraint qualification at
$x_{k}$ and (15) satisfies it at $s=0$, $x_{\star}$ is first-order critical.

We continue to follow the analysis of [3] and note that
$\Delta_{k}\geq\Delta_{\min}$ for all $k\in\mathds{N}$, where
$\Delta_{\min}:=\min(\Delta_{0},\,\gamma_{1}\Delta_{\textup{succ}})$. Aravkin
et al. [3] establish that $\xi_{\textup{cp}}(\cdot;x_{k},\nu_{k})$ is
increasing, and therefore,
$\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})\geq\xi_{\textup{cp}}(\Delta_{\min};x_{k},\nu_{k})$.
In addition, they show that
$\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}$ is an
appropriate criticality measure. Let $0<\epsilon<1$, and

$\displaystyle I(\epsilon)$
$\displaystyle:=\\{k\in\mathds{N}\mid\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}>\epsilon\\},$
$\displaystyle S(\epsilon)$ $\displaystyle:=\\{k\in
I(\epsilon)\mid\rho_{k}\geq\eta_{1}\\},$ $\displaystyle U(\epsilon)$
$\displaystyle:=\\{k\in I(\epsilon)\mid\rho_{k}<\eta_{1}\\},$

be the set of iterations, successful iterations, and unsuccessful iterations
until the criticality measure drops below $\epsilon$, respectively.

We now derive bounds on $|I(\epsilon)|$ using the analysis of [3] specialized
to Algorithm 5.1 under the assumption that $\\{D_{k}\\}$ is bounded. We make
the following assumption.

######  Assumption 1.

There exists $d_{\max}>0$ such that $\|d_{k}\|_{\infty}\leq d_{\max}$ for all
$k\in\mathds{N}$.

In Algorithm 5.1, it is not difficult to ensure satisfaction of 1. For
instance, one may prescribe a value $d_{\max}>0$ and reset each $d_{k}$
componentwise to $\min(\max(d_{k},-d_{\max}),d_{\max})$. We point out that 1
holds if (23) in Proposition 2 is satisfied for all $k$, because, in this
case, $\|d_{k}\|_{\infty}\leq
d_{\max}:=(\kappa_{\textup{mdc}}^{-1}-1)\alpha^{-1}\Delta_{\min}^{-1}$.

Under 1, the choice of $\nu_{k}$ in Algorithm 5.1 guarantees that

$\nu_{k}\geq\nu_{\min}:=1/(d_{\max}+\alpha^{-1}\Delta_{\min}^{-1})>0.$

The following result establishes that $|S(\epsilon)|$ is $O(\epsilon^{-2})$.

######  Proposition 5 (3, Lemma $3.6$).

Let Step Assumption 5.1 and 1 be satisfied. Assume that Algorithm 5.1
generates infinitely many successful iterations and that there exists
$(f+h)_{\textup{low}}\in\mathds{R}$ such that
$(f+h)(x_{k})\geq(f+h)_{\textup{low}}$ for all $k\in\mathds{N}$. Let
$\epsilon\in(0,\,1)$. Then,

$|S(\epsilon)|\leq\frac{(f+h)(x_{0})-(f+h)_{\textup{low}}}{\eta_{1}\kappa_{\textup{mdc}}\nu_{\min}\epsilon^{2}}.$

The next result establishes that $|U(\epsilon)|$ is also $O(\epsilon^{-2})$.

######  Proposition 6 (3, Lemma $3.7$).

Under the assumptions of Proposition 5,

$|U(\epsilon)|\leq\log_{\gamma_{2}}(\Delta_{\min}/\Delta_{0})+|S(\epsilon)|\,|\log_{\gamma_{2}}(\gamma_{4})|.$

The direct consequence of Propositions 5 and 6 is that $|I(\epsilon)|$ is also
$O(\epsilon^{-2})$. Because $\epsilon\in(0,\,1)$ is arbitrary, we conclude
that
$\liminf\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}=0$.

###  5.2 A variant saving proximal operator computations

We present a variant of Algorithm 5.1 in which we set

$\nu_{k}:=1/(\|d_{k}\|_{\infty}+\alpha^{-1})$ (25)

in 5, we do not compute $s_{k,1}$ in 7, and we leave $\Delta_{k}$ unchanged in
8, thus saving a proximal operator computation at each iteration. Under 1,
this choice of $\nu_{k}$ leads to a new $\nu_{\min}$

$\nu_{k}\geq\nu_{\min}:=1/(d_{\max}+\alpha^{-1})>0.$ (26)

We use the criticality measure

$\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2},$ (27)

and modify the definition of $I(\epsilon)$, $S(\epsilon)$ and $U(\epsilon)$
accordingly. To establish the same convergence properties as those of
Algorithm 5.1, we need Step Assumption 5.1 to hold for our algorithmic
variant. The following assumption summarizes our requirements.

######  Step Assumption 5.2.

There exists $\kappa_{\textup{m}}>0$ and $\kappa_{\textup{mdc}}\in(0,1)$ such
that for all $k$,

$\displaystyle|f(x_{k}+s_{k})+h(x_{k}+s_{k})-(\varphi(s_{k};x_{k},D_{k})+\psi(s_{k};x_{k}))|$
$\displaystyle\leq\kappa_{\textup{m}}\|s_{k}\|_{2}^{2},$ (28a)
$\displaystyle\xi(\Delta_{k};x_{k},D_{k})$
$\displaystyle\geq\kappa_{\textup{mdc}}\xi_{\textup{cp}}(\Delta_{k};x_{k};\nu_{k}),$
(28b)

where $\nu_{k}=1/(\|d_{k}\|_{\infty}+\alpha^{-1})$ in 5 and $s_{k}$ is
computed without changing $\Delta_{k}$ in 8 of Algorithm 5.1.

Even though Step Assumption 5.2 involves
$\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})$, we do not need to compute it in
practice. Conveniently, if $\\{D_{k}\\}$ is bounded, $\nabla f$ is Lipschitz
continuous, and we select $\psi(s;x_{k})=h(x_{k}+s)$, (28a) is still
satisfied. The following proposition, analogous to Proposition 2, shows that
(29) is a sufficient condition for (28b) to hold.

######  Proposition 7.

If

$\|d_{k}\|_{\infty}\leq(\kappa_{\textup{mdc}}^{-1}-1)\alpha^{-1},$ (29)

then (28b) is satisfied for Algorithm 5.1 with
$\nu_{k}=1/(\|d_{k}\|_{\infty}+\alpha^{-1})$ in 5, $s_{k,1}$ is not computed
in 7, and $\Delta_{k}$ is unchanged in 8.

###### Proof.

The proof is identical to that of Proposition 2 if we use (25). ∎

Whenever (28b) holds, we have

$\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},\nu_{k})^{1/2}\geq\kappa_{\textup{mdc}}^{1/2}\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}\geq\kappa_{\textup{mdc}}^{1/2}\nu_{k}^{-1}\|s_{k,1}\|,$
(30)

which ensures that
$\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},\nu_{k})^{1/2}\rightarrow 0$ implies
$\|s_{k,1}\|\rightarrow 0$, and that
$\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},\nu_{k})^{1/2}$ is an appropriate
criticality measure. We emphasize that there is no need to compute $s_{k,1}$
and $\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})$; (30) is only used as a
theoretical justification for the choice of our criticality measure.

We now establish convergence properties similar to those of Section 5.1, and
follow the analysis of Cartis et al. [13, Chapter $2.3$].

######  Proposition 8.

Let 1 and Step Assumption 5.2 be satisfied. Let $\nu_{\min}$ be as defined in
(26). If

$\Delta_{k}\leq\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2}$
(31)

and $x_{k}$ is not first-order stationary, iteration $k$ is very successful
and $\Delta_{k+1}>\Delta_{k}$.

###### Proof.

Using Step Assumption 5.2 and $\|s_{k}\|\leq\Delta_{k}$, we have

$\displaystyle|\rho_{k}-1|$
$\displaystyle=\left|\frac{f(x_{k}+s_{k})+h(x_{k}+s_{k})-m(s_{k};x_{k},D_{k})}{m(0;x_{k},D_{k})-m(s_{k};x_{k},D_{k})}\right|$
$\displaystyle\leq\frac{\kappa_{\textup{m}}\|s_{k}\|^{2}}{\xi(\Delta_{k};x_{k},D_{k})}$
$\displaystyle\leq\frac{\kappa_{\textup{m}}\Delta_{k}^{2}}{\xi(\Delta_{k};x_{k},D_{k})}.$

1 and (26) lead to

$|\rho_{k}-1|\leq\frac{\kappa_{\textup{m}}\nu_{\min}^{-1}\Delta_{k}^{2}}{\nu_{k}^{-1}\xi(\Delta_{k};x_{k},D_{k})}.$

Whenever (31) holds, we have $\rho_{k}\geq\eta_{2}$, implying that iteration
$k$ is very successful and $\Delta_{k+1}\geq\Delta_{k}$. ∎

The following lemma is inspired by Conn et al. [14, Theorem $6.4.3$] and
Cartis et al. [13, Lemma $2.3.4$].

######  Lemma 1.

Let Step Assumption 5.2 be satisfied. Then, for all $k\geq 0$,

$\Delta_{k}>\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(1,\Delta_{0}\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\min_{i\in[0,k]}\nu_{i}^{-1/2}\xi(\Delta_{i};x_{i},D_{i})^{1/2}.$
(32)

###### Proof.

We proceed as in Cartis et al. [13, Lemma $2.3.4$]. The bound certainly holds
for $k=0$, because

$\Delta_{0}>\gamma_{1}\sqrt{1-\eta_{2}}\Delta_{0}.$

Now, we proceed by contradiction and assume that $k\geq 1$ is the first
iteration such that (32) is not verified.

$\displaystyle\Delta_{k}$
$\displaystyle\leq\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(1,\Delta_{0}\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\min_{i\in[0,k]}\nu_{i}^{-1/2}\xi(\Delta_{i};x_{i},D_{i})^{1/2}$
$\displaystyle\leq\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\nu_{k-1}^{-1/2}\xi(\Delta_{k-1};x_{k-1},D_{k-1})^{1/2}.$

We have $\gamma_{1}\Delta_{k-1}\leq\Delta_{k}$ because of the updating rules
of $\Delta_{k}$ in Algorithm 5.1, which implies that

$\Delta_{k-1}\leq\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\nu_{k-1}^{-1/2}\xi(\Delta_{k-1};x_{k-1},D_{k-1})^{1/2}.$

Using Proposition 8, this results in iteration $k-1$ being very successful and
$\Delta_{k}>\Delta_{k-1}$. Therefore,

$\displaystyle\Delta_{k-1}<\Delta_{k}$
$\displaystyle\leq\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(1,\Delta_{0}\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\min_{i\in[0,k]}\nu_{i}^{-1/2}\xi(\Delta_{i};x_{i},D_{i})^{1/2}$
$\displaystyle\leq\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(1,\Delta_{0}\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\min_{i\in[0,k-1]}\nu_{i}^{-1/2}\xi(\Delta_{i};x_{i},D_{i})^{1/2},$

which is a contradiction with iteration $k$ being the first to violate (32). ∎

Now, we show analogous properties to Proposition 5 and Proposition 6 for our
variant.

######  Lemma 2.

Let Step Assumption 5.2 and 1 be satisfied. Assume that Algorithm 5.1 with
$\nu_{k}$ as in (25) in 5, without 7, and with $\Delta_{k}$ unchanged in 8
generates infinitely many successful iterations and that there exists
$(f+h)_{\textup{low}}\in\mathds{R}$ such that
$(f+h)(x_{k})\geq(f+h)_{\textup{low}}$ for all $k\in\mathds{N}$. Let
$\epsilon\in(0,\,1)$. Then,

$|S(\epsilon)|\leq\frac{(f+h)(x_{0})-(f+h)_{\textup{low}}}{\eta_{1}\nu_{\min}\epsilon^{2}}.$

###### Proof.

For $k\in S(\epsilon)$, we have

$\displaystyle f(x_{k})+h(x_{k})-f(x_{k}+s_{k})-h(x_{k}+s_{k})$
$\displaystyle\geq\eta_{1}(m(0;x_{k},D_{k})-m(s_{k};x_{k},D_{k}))$
$\displaystyle=\eta_{1}\xi(\Delta_{k};x_{k},D_{k})$
$\displaystyle\geq\eta_{1}\nu_{k}\epsilon^{2}$
$\displaystyle\geq\eta_{1}\nu_{\min}\epsilon^{2}.$

As $(f+h)(x_{k})\geq(f+h)_{\textup{low}}$, we sum the above inequalities for
$k\in S(\epsilon)$ and get

$(f+h)(x_{0})-(f+h)_{\textup{low}}\geq\sum_{k\in
S(\epsilon)}(f+h)(x_{k})-(f+h)(x_{k+1})\geq|S(\epsilon)|\eta_{1}\epsilon^{2}\nu_{\min}.$
∎

∎

######  Lemma 3.

Under the assumptions of Lemma 2,

$|U(\epsilon)|=O(\epsilon^{-2}).$

###### Proof.

Let $k<k(\epsilon)$. Because

$\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2}>\epsilon,$

and

$\nu_{k(\epsilon)}^{-1/2}\xi(\Delta_{k(\epsilon)};x_{k(\epsilon)},D_{k(\epsilon)})^{1/2}\leq\epsilon,$

we have with Lemma 1

$\displaystyle\Delta_{k(\epsilon)-1}$
$\displaystyle>\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(1,\Delta_{0}\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\min_{i\in[0,k]}\nu_{i}^{-1}\xi(\Delta_{i};x_{i},D_{i})$
$\displaystyle>\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(1,\Delta_{0}\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\epsilon.$

For each successful iteration, $\Delta_{k+1}\leq\gamma_{4}\Delta_{k}$, and for
each unsuccessful iteration, $\Delta_{k+1}\leq\gamma_{2}\Delta_{k}$, which
implies

$\Delta_{k(\epsilon)-1}\leq\Delta_{0}\gamma_{2}^{|U(\epsilon)|}\gamma_{4}^{|S(\epsilon)|}.$

By taking the natural logarithm of the above inequality, we have

$\displaystyle|U(\epsilon)|\log(\gamma_{2})+|S(\epsilon)|\log(\gamma_{4})$
$\displaystyle\geq\log(\Delta_{k(\epsilon)-1}/\Delta_{0})$
$\displaystyle\geq\log\left(\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(\Delta_{0}^{-1},\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\epsilon\right),$

so that, using the previous inequalities and Lemma 2

$\displaystyle|U(\epsilon)|$
$\displaystyle\leq\log_{\gamma_{2}}\left(\gamma_{1}\sqrt{\frac{\nu_{\min}(1-\eta_{2})}{\kappa_{\textup{m}}}}\min\left(\Delta_{0}^{-1},\frac{\sqrt{\kappa_{\textup{m}}\nu_{\min}^{-1}}}{\nu_{0}^{-1/2}\xi(\Delta_{0};x_{0},D_{0})^{1/2}}\right)\epsilon\right)+|S(\epsilon)||\log_{\gamma_{2}}(\gamma_{4})|$
$\displaystyle\leq O(|\log_{\gamma_{2}}(\epsilon)|)+O(\epsilon^{-2})$
$\displaystyle=O(\epsilon^{-2}).$ ∎

∎

Lemma 2 and Lemma 3 also indicate that $I(\epsilon)$ is $O(\epsilon^{-2})$,
and $\liminf\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2}=0$. With (28b) in
Step Assumption 5.2, we also have
$\liminf\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}=0$.

Finally, we emphasize the changes made in this section to Algorithm 5.1:

  * •

in 5, $\nu_{k}=1/(\|d_{k}\|_{\infty}+\alpha^{-1}\Delta_{k}^{-1})$ became
$\nu_{k}=1/(\|d_{k}\|_{\infty}+\alpha^{-1})$,

  * •

7 of Algorithm 5.1 was removed,

  * •

$\Delta_{k}$ was left unchanged in 8 of Algorithm 5.1,

  * •

we used $\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2}$ as a criticality
measure, instead of
$\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})$, however, when
(28a) holds, (30) indicates that if the new criticality measure is small, then
the criticality measure of Section 5.1 is also small,

  * •

in Proposition 8, there is no $\Delta_{\min}>0$ independent of $k$ as in
Proposition 3,

  * •

under the assumptions of Lemma 2, which are similar to those of Proposition 5,
the complexity bound is still in $O(\epsilon^{-2})$, and
$\liminf\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2}=\liminf\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}=0$.

##  6 Implementation and numerical experiments

Our Julia implementation of Algorithm 5.1 is available from our
RegularizedOptimization package [5] under the name TRDH. The latter can be
used to solve (1) directly, or as subproblem solver in TR [3, Algorithm
$3.1$], instead of R2 [3, Algorithm $6.1$]. Below, we use the notation TR-R2
and TR-TRDH to denote the application of TR to solve (1) with R2 or TRDH as
subproblem solver, respectively.

In our experiments, TR uses either an LSR1 or an LBFGS quasi-Newton Hessian
approximation with memory $5$, as implemented in the LinearOperators package
[25]. The same package implements diagonal quasi-Newton operators in TRDH
using the spectral update, the PSB update (12), and the update of Andrei [1]
(14) with the modification that we scaled the weak secant equations as

$\tilde{s}_{k-1}^{T}B\tilde{s}_{k-1}=\tilde{s}_{k-1}^{T}\tilde{y}_{k-1}$ (34)

in (11) and (13), where $\tilde{s}_{k-1}:=s_{k-1}/\|s_{k-1}\|_{2}$, and
$\tilde{y}_{k-1}:=y_{k-1}/\|s_{k-1}\|_{2}$, in order to alleviate numerical
issues as $s_{k-1}$ approaches zero.

The indefinite proximal operators of Example 1 and Example 2 are implemented
as part of the ShiftedProximalOperators package [7].

When using TRDH as the main solver, we initialize $D_{0}:=\nu_{0}^{-1}I$ for
$\nu_{0}>0$ given below. When using TR-TRDH with a spectral diagonal quasi-
Newton approximation, denoted TR-TRDH-Spec, we set the initial diagonal
Hessian approximation in TRDH at iteration $k$ of TR to
$D_{k,0}:=\nu_{k}^{-1}I$ (as we would initialize R2 in TR-R2). When using TR-
TRDH with the PSB or the Andrei quasi-Newton approximations, denoted TR-TRDH-
PSB and TR-TRDH-Andrei, respectively, we set
$D_{k,0}:=\mathop{\text{diag}}(B_{k})$, where $B_{k}$ is the quasi-Newton
Hessian approximation at iteration $k$ of TR.

We set $\psi(s;x_{k}):=h(x_{k}+s)$. We initialize $\nu_{0}=1$ for R2 and TRDH
used by themselves. The stopping criteria that we used for TR, TRDH and R2 (as
subproblem solvers or main solvers) are based on
$\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}$. We set $\Delta_{0}=1$ for
TR and for TRDH used as main solver. For TR-TRDH, at iteration $k$ of TR, the
initial value of the TRDH trust-region radius is
$\Delta_{k,0}=\min(\Delta_{k},\beta\|s_{k,1}\|)$ / 10, where $\Delta_{k}$ is
the TR trust-region radius at iteration $k$, and $s_{k,1}$ is the first step
of the $k$-th TR subproblem. In other words, the initial TRDH trust-region
radius to solve the $k$-th TR subproblem is a tenth of the trust-region radius
of this $k$-th subproblem.

For all solvers, the outer iterations terminate as soon as

$\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}<\epsilon_{a}+\epsilon_{r}\nu_{0}^{-1/2}\xi_{\textup{cp}}(\Delta_{0};x_{0},\nu_{0})^{1/2},$
(35)

where $\epsilon_{a}>0$ and $\epsilon_{r}>0$ are an absolute and a relative
tolerance. A round of inner iterations in TR terminates as soon as the
stationarity measure of the inner solver satisfies (35), with
$\tilde{\epsilon}_{a}=10^{-5}$ for the first TR iteration, otherwise
$\tilde{\epsilon}_{a}=\max(\epsilon_{a,i},\min(10^{-2},\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}))$,
and $\tilde{\epsilon}_{r}=\epsilon_{r,i}$, where $\epsilon_{a,i}>0$ and
$\epsilon_{r,i}>0$ are some absolute and relative inner tolerances, and
$s_{k,1}$ is the first iterate of the solution of the trust-region subproblem.
In the experiments below, we use $\epsilon_{a,i}=10^{-3}$ and
$\epsilon_{r,i}=10^{-6}$, except in Section 6.1, where we use
$\epsilon_{a,i}=10^{-5}$.

Additionally, we test the variant presented in Section 5.2, which is denoted
“iTRDH” (_indefinite trust-region with diagonal Hessian approximations_) in
our results. When in use, the stopping criterion is based on
$\nu_{k}^{-1/2}\xi(\Delta_{k};x_{k},D_{k})^{1/2}$, instead of
$\nu_{k}^{-1/2}\xi_{\textup{cp}}(\Delta_{k};x_{k},\nu_{k})^{1/2}$.

In our results, we report

  * •

the final $f(x)$;

  * •

the final $h(x)/\lambda$;

  * •

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

  * •

$\|x-x_{\star}\|_{2}$, where $x_{\star}$ is the exact solution, if it is
available;

  * •

the number of smooth objective evaluations $\\#~{}f$;

  * •

the number of gradient evaluations $\\#~{}\nabla f$;

  * •

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

  * •

the elapsed time $t$ in seconds.

Because our implementations are not yet perfectly optimized in terms of memory
allocations, we neglect the elapsed time in our interpretations, and only
report it in the tables as an indicator for the reader.

In our test cases, available from the RegularizedProblems package [6], the
computational cost of evaluating the gradient is significantly higher than a
proximal evaluation or an objective evaluation. In all cases, except in case
of failure, all solvers find similar final solutions, and we only show one for
illustration.

###  6.1 Basis pursuit denoise (BPDN)

Our first test case is the basis pursuit denoise (BPDN) problem [28, 16]. The
stopping tolerances $\epsilon_{a}$ and $\epsilon_{r}$ are set to $10^{-5}$. In
this subsection only, when using TR, we set $\epsilon_{a,i}=10^{-5}$ in order
to have accurate subproblem solves, which leads to performing fewer gradient
evaluations without sacrificing too many proximal evaluations. 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 $10$ of its components that are set to $\pm 1$. We solve

$\mathop{\textup{minimize}}_{x}\tfrac{1}{2}\|Ax-b\|_{2}^{2}+h(x),$ (36)

where $h(x)=\lambda\|x\|_{0}$. As in [3], we use
$\lambda=0.1\|A^{T}b\|_{\infty}$.

Table 1:  BPDN (36) statistics with $h=\lambda\|\cdot\|_{0}$. All variants of TR use an LSR1 Hessian approximation, and are given a maximum of $100$ inner iterations. The optimal objective value is $f(x_{\star})=9.90e-03$. solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\|x-x_{\star}\|_{2}$ | $\\#f$ | $\\#\nabla f$ | $\\#\mathop{\textup{prox}}{}$ |  $t$ ($s$)
---|---|---|---|---|---|---|---|---
R2 |  $9.44$e$-03$ | $10$ |  $4.4$e$-03$ |  $4.7$e$-02$ | $30$ | $31$ | $30$ |  $8.0$e$-03$
TRDH-Spec |  $9.44$e$-03$ | $10$ |  $4.0$e$-03$ |  $4.7$e$-02$ | $9$ | $9$ | $17$ |  $1.7$e$-02$
iTRDH-Spec |  $9.44$e$-03$ | $10$ |  $3.4$e$-03$ |  $4.7$e$-02$ | $10$ | $9$ | $9$ |  $4.0$e$-03$
TRDH-PSB |  $9.44$e$-03$ | $10$ |  $3.7$e$-03$ |  $4.7$e$-02$ | $18$ | $15$ | $35$ |  $1.1$e$-02$
iTRDH-PSB |  $9.44$e$-03$ | $10$ |  $2.6$e$-03$ |  $4.7$e$-02$ | $20$ | $16$ | $19$ |  $5.0$e$-03$
TRDH-Andrei |  $9.44$e$-03$ | $10$ |  $4.5$e$-03$ |  $4.7$e$-02$ | $48$ | $29$ | $95$ |  $1.0$e$-02$
iTRDH-Andrei |  $9.44$e$-03$ | $10$ |  $3.7$e$-03$ |  $4.7$e$-02$ | $59$ | $35$ | $58$ |  $3.1$e$-02$
TR-R2 |  $9.44$e$-03$ | $10$ |  $1.7$e$-05$ |  $4.7$e$-02$ | $23$ | $23$ | $40$ |  $1.8$e$-02$
TR-TRDH-PSB |  $9.44$e$-03$ | $10$ |  $1.6$e$-05$ |  $4.7$e$-02$ | $21$ | $21$ | $58$ |  $6.9$e$-02$
TR-iTRDH-PSB |  $9.44$e$-03$ | $10$ |  $1.6$e$-05$ |  $4.7$e$-02$ | $21$ | $21$ | $39$ |  $2.5$e$-02$
TR-TRDH-Andrei |  $9.44$e$-03$ | $10$ |  $1.6$e$-05$ |  $4.7$e$-02$ | $20$ | $20$ | $127$ |  $2.5$e$-02$
TR-iTRDH-Andrei |  $9.44$e$-03$ | $10$ |  $1.4$e$-05$ |  $4.7$e$-02$ | $20$ | $20$ | $100$ |  $3.0$e$-02$
TR-TRDH-Spec |  $9.44$e$-03$ | $10$ |  $1.9$e$-05$ |  $4.7$e$-02$ | $21$ | $21$ | $46$ |  $1.8$e$-02$
TR-iTRDH-Spec |  $9.44$e$-03$ | $10$ |  $1.9$e$-05$ |  $4.7$e$-02$ | $21$ | $21$ | $33$ |  $1.8$e$-02$

Figure 3:  Solution of (36) (left), and error (right) with TRDH-Spec.

Figure 3 shows the solution of (36) with TRDH-Spec. Table 1 reports the
statistics for the different solvers and shows that TRDH-Spec and TRDH-PSB
perform fewer objective and gradient evaluations than R2. TRDH-Andrei performs
worse that R2 on this problem. As expected, the “iTRDH” variants require fewer
proximal operator evaluations. They result in similar numbers of objective and
gradient evaluations in the cases of TRDH-Spec and TRDH-PSB, but require more
evaluations in the case of TRDH-Andrei.

All TR-TRDH and TR-iTRDH solvers perform fewer objective and gradient
evaluations than TR-R2.

We now solve the constrained variant

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

where each element in $x_{\star}$ is either $0$ or $1$.

Table 2:  Constrained BPDN (37) statistics with $h=\lambda\|\cdot\|_{0}$. All variants of TR use an LSR1 Hessian approximation, and are given a maximum of $100$ inner iterations. The optimal objective value is $f(x_{\star})=8.98e-03$. solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\|x-x_{\star}\|_{2}$ | $\\#f$ | $\\#\nabla f$ | $\\#\mathop{\textup{prox}}{}$ |  $t$ ($s$)
---|---|---|---|---|---|---|---|---
R2 |  $8.70$e$-03$ | $10$ |  $3.9$e$-03$ |  $3.6$e$-02$ | $27$ | $28$ | $27$ |  $1.9$e$-02$
TRDH-Spec |  $8.70$e$-03$ | $10$ |  $3.6$e$-03$ |  $3.6$e$-02$ | $8$ | $8$ | $15$ |  $4.0$e$-03$
iTRDH-Spec |  $8.70$e$-03$ | $10$ |  $3.0$e$-03$ |  $3.6$e$-02$ | $9$ | $8$ | $8$ |  $2.0$e$-03$
TRDH-PSB |  $8.70$e$-03$ | $10$ |  $3.6$e$-03$ |  $3.6$e$-02$ | $11$ | $11$ | $21$ |  $3.0$e$-03$
iTRDH-PSB |  $8.70$e$-03$ | $10$ |  $3.7$e$-03$ |  $3.6$e$-02$ | $12$ | $11$ | $11$ |  $3.0$e$-03$
TRDH-Andrei |  $8.70$e$-03$ | $10$ |  $4.3$e$-03$ |  $3.6$e$-02$ | $48$ | $29$ | $95$ |  $1.6$e$-02$
iTRDH-Andrei |  $8.70$e$-03$ | $10$ |  $4.0$e$-03$ |  $3.6$e$-02$ | $59$ | $35$ | $58$ |  $2.0$e$-02$
TR-R2 |  $8.70$e$-03$ | $10$ |  $1.9$e$-05$ |  $3.6$e$-02$ | $20$ | $20$ | $32$ |  $1.2$e$-02$
TR-TRDH-PSB |  $8.70$e$-03$ | $10$ |  $2.1$e$-05$ |  $3.6$e$-02$ | $19$ | $19$ | $50$ |  $1.8$e$-02$
TR-iTRDH-PSB |  $8.70$e$-03$ | $10$ |  $2.1$e$-05$ |  $3.6$e$-02$ | $19$ | $19$ | $34$ |  $2.5$e$-02$
TR-TRDH-Andrei |  $8.70$e$-03$ | $10$ |  $1.9$e$-05$ |  $3.6$e$-02$ | $20$ | $20$ | $63$ |  $1.8$e$-02$
TR-iTRDH-Andrei |  $8.70$e$-03$ | $10$ |  $1.5$e$-05$ |  $3.6$e$-02$ | $20$ | $20$ | $51$ |  $7.0$e$-03$
TR-TRDH-Spec |  $8.70$e$-03$ | $10$ |  $2.0$e$-05$ |  $3.6$e$-02$ | $19$ | $19$ | $28$ |  $1.0$e$-02$
TR-iTRDH-Spec |  $8.70$e$-03$ | $10$ |  $2.0$e$-05$ |  $3.6$e$-02$ | $19$ | $19$ | $23$ |  $1.2$e$-02$

Figure 4:  Solution of (37) (left), and error (right) with TR-TRDH-Spec.

Figure 4 shows the solution of (37) with TRDH-Spec, and Table 2 reports the
statistics for the different solvers. We observe that TRDH-Spec, iTRDH-Spec,
TRDH-PSB and iTRDH-PSB used as main solvers outperform R2 in terms of number
of objective and gradient evaluations. TR-TRDH and TR-iTRDH perform similarly
to TR-R2.

###  6.2 Sparse nonnegative matrix factorization (NNMF)

Our next test case is a variant of the NNMF problem of Kim and Park [20]. Let
$A\in\mathds{R}^{m\times n}$ have nonnegative entries, where each row
represents a feature and each column represents an observation. We wish to
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. The problem is stated as

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

where $h(H)=\lambda\|\mathop{\textup{vec}}{H}\|_{0}$ and
$\mathop{\textup{vec}}{}$ stacks the columns of a matrix to form a vector. In
our experiments, each observation is generated using a mixture of Gaussians.
Negative elements in the matrix $A$ generated are reset to zero.

Table 3:  NNMF (38) statistics with $h=\lambda\|\cdot\|_{0}$. All variants of TR use an LSR1 Hessian approximation, and are given a maximum of $100$ inner iterations. solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\\#f$ | $\\#\nabla f$ | $\\#\mathop{\textup{prox}}{}$ |  $t$ ($s$)
---|---|---|---|---|---|---|---
R2 |  $2.84$e$+03$ | $0$ |  $2.2$e$-04$ | $2$ | $2$ | $2$ |  $0.0$e$+00$
TRDH-Spec |  $1.25$e$+02$ | $50$ |  $7.7$e$-02$ | $47$ | $29$ | $93$ |  $5.0$e$-03$
iTRDH-Spec |  $1.25$e$+02$ | $50$ |  $6.5$e$-02$ | $48$ | $29$ | $47$ |  $4.0$e$-03$
TRDH-PSB |  $1.70$e$+02$ | $62$ |  $5.7$e$+00$ | $501$ | $352$ | $1000$ |  $5.6$e$-02$
iTRDH-PSB |  $1.70$e$+02$ | $62$ |  $8.4$e$+00$ | $501$ | $352$ | $500$ |  $5.5$e$-02$
TRDH-Andrei |  $2.73$e$+02$ | $54$ |  $5.9$e$+00$ | $501$ | $283$ | $1000$ |  $5.9$e$-02$
iTRDH-Andrei |  $2.73$e$+02$ | $54$ |  $7.0$e$+00$ | $501$ | $283$ | $500$ |  $5.6$e$-02$
TR-R2 |  $1.25$e$+02$ | $50$ |  $5.0$e$-03$ | $212$ | $115$ | $4983$ |  $2.4$e$-01$
TR-TRDH-PSB |  $1.25$e$+02$ | $50$ |  $7.2$e$-03$ | $197$ | $104$ | $9826$ |  $7.3$e$-01$
TR-iTRDH-PSB |  $1.25$e$+02$ | $50$ |  $6.4$e$-03$ | $129$ | $71$ | $3759$ |  $4.8$e$-01$
TR-TRDH-Andrei |  $1.25$e$+02$ | $50$ |  $2.9$e$-03$ | $170$ | $91$ | $8746$ |  $6.1$e$-01$
TR-iTRDH-Andrei |  $1.25$e$+02$ | $50$ |  $6.8$e$-03$ | $152$ | $83$ | $4809$ |  $6.1$e$-01$
TR-TRDH-Spec |  $1.25$e$+02$ | $50$ |  $4.8$e$-03$ | $126$ | $69$ | $4793$ |  $3.6$e$-01$
TR-iTRDH-Spec |  $1.25$e$+02$ | $50$ |  $3.6$e$-03$ | $162$ | $82$ | $2503$ |  $3.5$e$-01$

Figure 5:  Solution of (38) with TR-TRDH-PSB.

We set $m=100$, $n=50$, $k=5$, $\lambda=10^{-1}$. The stopping tolerances
$\epsilon_{a}$ and $\epsilon_{r}$ are set to $10^{-5}$. Figure 5 shows the
solution of (38) with TR-TRDH-PSB. The statistics are reported in Table 3. We
observe that R2 is trapped in a spurious stationary point, and that all TRDH
and iTRDH solvers except TRDH-Spec and iTRDH-Spec reach their maximum number
of iterations. All solvers using TR perform well with fewer objective and
gradient evaluations than TR-R2. The number of proximal operator calls is
lower for all TR-iTRDH variants.

###  6.3 Nonlinear support vector machine (SVM)

We now consider the nonlinear SVM described by [2] to classify digits of the
MNIST dataset as either $1$ or $7$—the other digits are removed. Let $m$ be
the number of images and $n$ the vectorized image size. The problem reads

$\mathop{\textup{minimize}}_{x}\tfrac{1}{2}\|\textbf{1}-\tanh(b\odot(A^{T}x))\|^{2}+h(x),$
(39)

where $A\in\mathds{R}^{m\times n}$ is the data matrix, $b$ is the vector of
labels with values $\pm 1$ for all its components, $\odot$ denotes the
elementwise product between two vectors, and $h(x)=\lambda\|x\|_{1}$. We use
$n=784$, $m=13007$ in the training set, $m=2163$ in the testing set,
$\lambda=0.1$ and we initialize the problem at $x=e$, the vector of ones. We
set the absolute and relative stopping tolerances $\epsilon_{a}$ and
$\epsilon_{r}$ to $10^{-4}$.

Table 4:  SVM (39) statistics with $h=\lambda\|\cdot\|_{1}$. All variants of TR use an LBFGS Hessian approximation, and are given a maximum of $100$ inner iterations. The train and test accuracies are the percentage of images correctly classified, computed by counting the number of elements of the residual $\textbf{1}-\tanh(b\odot(A^{T}x))$ that are lower than $1$ for the train and test problems respectively.  solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | (Train, Test) | $\\#f$ | $\\#\nabla f$ | $\\#\mathop{\textup{prox}}{}$ |  $t$ ($s$)
---|---|---|---|---|---|---|---|---
R2 |  $2.18$e$+02$ |  $2.4$e$+03$ |  $1.2$e$-01$ | $(99.3,98.8)$ | $265$ | $199$ | $265$ |  $3.7$e$+00$
TRDH-Spec |  $2.18$e$+02$ |  $2.4$e$+03$ |  $5.9$e$-02$ | $(99.3,98.9)$ | $306$ | $191$ | $611$ |  $4.4$e$+00$
iTRDH-Spec |  $2.18$e$+02$ |  $2.4$e$+03$ |  $5.9$e$-02$ | $(99.3,98.9)$ | $307$ | $191$ | $306$ |  $3.9$e$+00$
TRDH-PSB |  $2.38$e$+02$ |  $2.4$e$+03$ |  $4.5$e$-01$ | $(99.1,98.8)$ | $1001$ | $686$ | $2000$ |  $1.3$e$+01$
iTRDH-PSB |  $2.38$e$+02$ |  $2.4$e$+03$ |  $4.3$e$-01$ | $(99.1,98.8)$ | $1001$ | $686$ | $1000$ |  $1.3$e$+01$
TRDH-Andrei |  $3.00$e$+02$ |  $2.9$e$+03$ |  $1.2$e$+00$ | $(99.1,99.0)$ | $1001$ | $418$ | $2000$ |  $9.1$e$+00$
iTRDH-Andrei |  $3.00$e$+02$ |  $2.9$e$+03$ |  $2.7$e$+00$ | $(99.1,99.0)$ | $1001$ | $418$ | $1000$ |  $8.8$e$+00$
TR-R2 |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.7$e$-03$ | $(99.3,98.8)$ | $361$ | $361$ | $1512$ |  $5.3$e$+00$
TR-TRDH-PSB |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.7$e$-03$ | $(99.3,98.8)$ | $316$ | $316$ | $3882$ |  $6.5$e$+00$
TR-iTRDH-PSB |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.7$e$-03$ | $(99.3,98.8)$ | $377$ | $377$ | $2242$ |  $7.8$e$+00$
TR-TRDH-Andrei |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.7$e$-03$ | $(99.3,98.9)$ | $396$ | $396$ | $27155$ |  $1.4$e$+01$
TR-iTRDH-Andrei |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.5$e$-03$ | $(99.3,98.9)$ | $283$ | $283$ | $11418$ |  $1.2$e$+01$
TR-TRDH-Spec |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.7$e$-03$ | $(99.3,98.8)$ | $333$ | $333$ | $3390$ |  $7.5$e$+00$
TR-iTRDH-Spec |  $2.18$e$+02$ |  $2.4$e$+03$ |  $4.7$e$-03$ | $(99.3,98.8)$ | $472$ | $472$ | $1791$ |  $9.3$e$+00$

Figure 6:  Solution map of (39) with TR-TRDH-Spec (left), and sample 1’s
(right, top row) and 7’s (right, bottom row) rotated digits from the MNIST
dataset.

The solution map of (39) representing the most important pixels to classify
the images between $1$ and $7$ is shown in Figure 6 for the solver TR-TRDH-
Spec, along with a few sample digits. We observe that the mid-height pixels
are selected as the most important ones, which is consistent with the fact
that the images show the digits sideways. The statistics are reported in Table
4, and show that TRDH-PSB, iTRDH-PSB, TRDH-Andrei and iTRDH-Andrei as main
solvers exceed the maximum number of iterations. R2, TRDH-Spec and iTRDH-Spec
are the most efficient and perform similar number of objective and gradient
evaluations, with R2 terminating with a slightly higher criticality measure.
TR-TRDH-PSB, TR-iTRDH-Andrei and TR-TRDH-Spec perform fewer objective and
gradient evaluations, but more proximal operator evaluations than TR-R2.

###  6.4 FitzHugh-Nagumo inverse problem

Let $v(x)=(v_{1}(x),\ldots,v_{n+1}(x))$ and
$w(x)=(w_{1}(x),\ldots,w_{n+1}(x))$ be sampled values of $V(t;x)$ and $W(t;x)$
for $x\in\mathds{R}^{5}$ satisfying the FitzHugh [18] and Nagumo et al. [23]
model for neuron activation

$\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}).$
(40)

The samples are generated using a discretization of the time interval
$t\in[0,20]$ with initial conditions $(V(0),W(0))=(2,0)$. We define a target
solution that corresponds to a simulation of the van der Pol [29] oscillator
by generating solutions $(\bar{v}(x),\bar{w}(x))$ of (40) with
$\bar{x}=(0,0.2,1,0,0)$, and we solve

$\mathop{\textup{minimize}}_{x}\tfrac{1}{2}\|(v(x)-\bar{v}(\bar{x}),w(x)-\bar{w}(\bar{x}))\|_{2}^{2}+h(x),$
(41)

where $h(x)=\lambda\|x\|_{0}$ with $n=100$. The stopping tolerances
$\epsilon_{a}$ and $\epsilon_{r}$ are set to $10^{-4}$. We report the
statistics of the solution of (41) with $\lambda=10$ in Table 5. TRDH as main
solver results in many objective and gradient evaluations compared to TR.
Moreover, TR-TRDH-PSB, TR-iTRDH-PSB, TR-TRDH-Andrei and TR-iTRDH-Andrei
perform fewer objective and gradient evaluations than TR-R2. TR-TRDH-Spec is
the least efficient.

Table 5:  FH (41) statistics with $h=\lambda\|\cdot\|_{0}$. All variants of TR use an LBFGS Hessian approximation, and are given a maximum of $200$ inner iterations. The optimal objective value is $f(x_{\star})=1.03e+00$. solver | $f(x)$ | $h(x)/\lambda$ | $\sqrt{\xi/\nu}$ | $\|x-x_{\star}\|_{2}$ | $\\#f$ | $\\#\nabla f$ | $\\#\mathop{\textup{prox}}{}$ |  $t$ ($s$)
---|---|---|---|---|---|---|---|---
TRDH-Spec |  $3.21$e$+01$ | $5$ |  $7.0$e$+00$ |  $9.8$e$-01$ | $501$ | $351$ | $1000$ |  $2.1$e$+00$
iTRDH-Spec |  $3.21$e$+01$ | $5$ |  $5.9$e$+00$ |  $9.8$e$-01$ | $501$ | $351$ | $500$ |  $2.5$e$+00$
TRDH-PSB |  $1.12$e$+00$ | $2$ |  $2.7$e$+00$ |  $1.3$e$-01$ | $501$ | $401$ | $1000$ |  $4.7$e$+00$
iTRDH-PSB |  $1.12$e$+00$ | $2$ |  $2.4$e$+00$ |  $1.3$e$-01$ | $501$ | $401$ | $500$ |  $4.7$e$+00$
TRDH-Andrei |  $1.15$e$+00$ | $2$ |  $1.6$e$+00$ |  $1.5$e$-01$ | $501$ | $372$ | $1000$ |  $4.2$e$+00$
iTRDH-Andrei |  $1.15$e$+00$ | $2$ |  $2.8$e$+00$ |  $1.5$e$-01$ | $501$ | $372$ | $500$ |  $5.3$e$+00$
TR-R2 |  $1.02$e$+00$ | $2$ |  $5.0$e$-03$ |  $6.3$e$-03$ | $327$ | $236$ | $36274$ |  $6.7$e$+00$
TR-TRDH-PSB |  $1.02$e$+00$ | $2$ |  $5.7$e$-03$ |  $6.1$e$-03$ | $183$ | $149$ | $56042$ |  $3.8$e$+00$
TR-iTRDH-PSB |  $1.02$e$+00$ | $2$ |  $6.3$e$-03$ |  $6.2$e$-03$ | $211$ | $179$ | $33452$ |  $3.7$e$+00$
TR-TRDH-Andrei |  $1.02$e$+00$ | $2$ |  $5.6$e$-03$ |  $6.5$e$-03$ | $185$ | $138$ | $49311$ |  $2.9$e$+00$
TR-iTRDH-Andrei |  $1.02$e$+00$ | $2$ |  $6.4$e$-03$ |  $6.1$e$-03$ | $168$ | $151$ | $23655$ |  $2.9$e$+00$
TR-TRDH-Spec |  $1.04$e$+00$ | $2$ |  $1.0$e$+00$ |  $7.3$e$-02$ | $501$ | $428$ | $109750$ |  $9.0$e$+00$
TR-iTRDH-Spec |  $1.02$e$+00$ | $2$ |  $4.6$e$-03$ |  $1.1$e$-02$ | $383$ | $337$ | $40354$ |  $6.1$e$+00$
Table 6:  FH (41) (left) and constrained FH (42) (right) solutions identified by the solvers tested. The unconstrained solution is $(0,0.2,1,0,0)$.  | $x_{1}$ | $x_{2}$ | $x_{3}$ | $x_{4}$ | $x_{5}$ | $x_{1}$ | $x_{2}$ | $x_{3}$ | $x_{4}$ | $x_{5}$
---|---|---|---|---|---|---|---|---|---|---
TRDH-Spec | 0.32 | 0.43 | 0.83 | 0.77 | 0.44 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
iTRDH-Spec | 0.32 | 0.43 | 0.83 | 0.77 | 0.44 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TRDH-PSB | 0.00 | 0.24 | 0.87 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
iTRDH-PSB | 0.00 | 0.24 | 0.87 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TRDH-Andrei | 0.00 | 0.25 | 0.86 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
iTRDH-Andrei | 0.00 | 0.25 | 0.86 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-R2 | 0.00 | 0.20 | 1.01 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-TRDH-PSB | 0.00 | 0.20 | 1.01 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-iTRDH-PSB | 0.00 | 0.20 | 1.01 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-TRDH-Andrei | 0.00 | 0.20 | 1.01 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-iTRDH-Andrei | 0.00 | 0.20 | 1.01 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-TRDH-Spec | 0.00 | 0.18 | 1.07 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00
TR-iTRDH-Spec | 0.00 | 0.20 | 1.01 | 0.00 | 0.00 | 0.00 | 0.50 | 0.54 | 0.00 | 0.00

Figure 7:  Solution of (41) (left) and of (42) (right) with TR-TRDH-PSB and
sampled values of $V$ and $W$.

The left half of Table 6 reports the final solution identified by each solver.
The solution of (41) with TR-TRDH-PSB and the sampled values of $V$ and $W$
are displayed in Figure 7. We do not show results with R2 because it
encountered numerical issues when solving the problem.

We also solve the constrained variant

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

and we keep all parameters the same, except for $\lambda$ that we set to $40$
to further enforce sparsity.

Table 7:  Constrained FH (42) statistics with $h=\lambda\|\cdot\|_{1}$. All variants of TR use an LBFGS Hessian approximation, and are given a maximum of $200$ inner iterations. solver | $f(x)$ | $h(x)/\lambda$ | $\xi$ | $\|x-x_{\star}\|_{2}$ | $\\#f$ | $\\#\nabla f$ | $\\#\mathop{\textup{prox}}{}$ |  $t$ ($s$)
---|---|---|---|---|---|---|---|---
TRDH-Spec |  $4.43$e$+00$ |  $1.0$e$+00$ |  $1.2$e$-02$ |  $5.5$e$-01$ | $327$ | $226$ | $653$ |  $2.6$e$+00$
iTRDH-Spec |  $4.43$e$+00$ |  $1.0$e$+00$ |  $1.1$e$-02$ |  $5.5$e$-01$ | $328$ | $226$ | $327$ |  $2.4$e$+00$
TRDH-PSB |  $4.43$e$+00$ |  $1.0$e$+00$ |  $8.2$e$-03$ |  $5.5$e$-01$ | $331$ | $237$ | $661$ |  $2.6$e$+00$
iTRDH-PSB |  $4.43$e$+00$ |  $1.0$e$+00$ |  $7.7$e$-03$ |  $5.5$e$-01$ | $332$ | $237$ | $331$ |  $2.2$e$+00$
TRDH-Andrei |  $4.43$e$+00$ |  $1.0$e$+00$ |  $7.1$e$-03$ |  $5.5$e$-01$ | $223$ | $155$ | $445$ |  $1.4$e$+00$
iTRDH-Andrei |  $4.43$e$+00$ |  $1.0$e$+00$ |  $7.1$e$-03$ |  $5.5$e$-01$ | $224$ | $155$ | $223$ |  $1.6$e$+00$
TR-R2 |  $4.43$e$+00$ |  $1.0$e$+00$ |  $4.9$e$-03$ |  $5.5$e$-01$ | $47$ | $31$ | $3057$ |  $4.4$e$-01$
TR-TRDH-PSB |  $4.43$e$+00$ |  $1.0$e$+00$ |  $2.7$e$-04$ |  $5.5$e$-01$ | $38$ | $29$ | $6725$ |  $3.1$e$-01$
TR-iTRDH-PSB |  $4.43$e$+00$ |  $1.0$e$+00$ |  $7.9$e$-05$ |  $5.5$e$-01$ | $36$ | $27$ | $3326$ |  $3.3$e$-01$
TR-TRDH-Andrei |  $4.43$e$+00$ |  $1.0$e$+00$ |  $8.7$e$-05$ |  $5.5$e$-01$ | $39$ | $30$ | $7438$ |  $3.8$e$-01$
TR-iTRDH-Andrei |  $4.43$e$+00$ |  $1.0$e$+00$ |  $1.6$e$-03$ |  $5.5$e$-01$ | $36$ | $27$ | $3117$ |  $3.4$e$-01$
TR-TRDH-Spec |  $4.43$e$+00$ |  $1.0$e$+00$ |  $4.1$e$-03$ |  $5.5$e$-01$ | $36$ | $27$ | $5467$ |  $3.9$e$-01$
TR-iTRDH-Spec |  $4.43$e$+00$ |  $1.0$e$+00$ |  $1.8$e$-03$ |  $5.5$e$-01$ | $38$ | $25$ | $3347$ |  $2.5$e$-01$

The right half of Table 6 reports the solution and Table 7 the statistics of
the solve of (42) with the tested solvers. Note that all solvers identify the
same solution. We observe that all TR-TRDH and TR-iTRDH solvers perform fewer
objective and gradient evaluations than TR-R2. However, TR-R2 has the lowest
number of proximal operator calls out of all solvers using TR. The right plot
of Figure 7 shows that the $W$ part of the solution of (42) with TR-TRDH-PSB
does not match the data as tightly as that of (41) in the left half of Figure
7. That is a consequence of enforcing $x_{2}\geq 0.5$, since the unconstrained
solution verifies $x_{2}=0.2$.

##  7 Discussion and future work

The $\ell_{0}$\- and $\ell_{1}$-norm regularizers are standard choices to
promote sparsity. It is possible to derive the $\mathop{\textup{iprox}}$ of
other useful separable regularizers, including $\ell_{p}$ pseudonorms to the
$p$-th power, i.e.,

$h(x)=\|x\|_{p}^{p}=\sum_{i=1}^{n}|x_{i}|^{p}\quad(0<p<1),$

and those are also useful to promote sparsity [12]. Much remains to be done,
however, including deriving the $\mathop{\textup{iprox}}$ of relevant non-
separable regularizers, and studying other diagonal Hessian approximations
than the ones considered above, including new diagonal quasi-Cauchy updates.

TRDH performs well on the problems tested and is promising, but we were
surprised to see the spectral gradient update often perform better than more
sophisticated diagonal quasi-Newton updates, especially when using TRDH as
main solver. Future research should seek to provide an explanation for that
observation. In several instances, the variant of TRDH denoted “iTRDH”
performs better than the basic version, which indicates that both algorithms
are relevant.

Finally, other solvers would likely benefit from using TRDH as a subproblem
solver, including the methods for least-squares $f$ of Aravkin et al. [2].

### Acknowledgements

The authors thank Joshua Wolff from École Normale Supérieure des Techniques
Avancées (ENSTA), Paris, for the work he conducted during his undergraduate
internship at GERAD in the summer of 2022 and that made this research
possible.

## References

  * Andrei [2019] N. Andrei.  A diagonal quasi-Newton updating method for unconstrained optimization.  _Numer. Algor._ , 81:575––590, 2019.
  * Aravkin et al. [2022a] A. Aravkin, R. Baraldi, and D. Orban.  A Levenberg-Marquardt method for nonsmooth regularized least squares.  Cahier du GERAD G-2023-58, GERAD, Montréal, QC, Canada, 2022a.
  * Aravkin et al. [2022b] 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, 2022b.
  * Baraldi and Kouri [2022] R. Baraldi and D. P. Kouri.  A proximal trust-region method for nonsmooth optimization with inexact function and gradient evaluations.  _Math. Program._ , 2022.
  * Baraldi and Orban [2022a] R. Baraldi and D. Orban.  RegularizedOptimization.jl: Algorithms for regularized optimization.  https://github.com/JuliaSmoothOptimizers/RegularizedOptimization.jl, February 2022a.
  * Baraldi and Orban [2022b] R. Baraldi and D. Orban.  RegularizedProblems.jl: Test cases for regularized optimization.  https://github.com/JuliaSmoothOptimizers/RegularizedProblems.jl, February 2022b.
  * Baraldi and Orban [2022c] R. Baraldi and D. Orban.  ShiftedProximalOperators.jl: Proximal operators for regularized optimization.  https://github.com/JuliaSmoothOptimizers/ShiftedProximalOperators.jl, February 2022c.
  * Beck [2017] A. Beck.  _First-Order Methods in Optimization_.  Number 25 in MOS-SIAM Series on Optimization. SIAM, Philadelphia, USA, 2017.
  * Becker and Fadili [2012] S. Becker and J. Fadili.  A quasi-Newton proximal splitting method.  In F. Pereira, C. Burges, L. Bottou, and K. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 25. Curran Associates, Inc., 2012.
  * Becker et al. [2019] S. Becker, J. Fadili, and P. Ochs.  On quasi-Newton forward-backward splitting: Proximal calculus and convergence.  _SIAM J. Optim._ , 29(4):2445–2481, 2019.
  * Birgin et al. [2014] E. G. Birgin, J. M. Martínez, and M. Raydan.  Spectral projected gradient methods: Review and perspectives.  _Journal of Statistical Software_ , 60(3):1–21, 2014.
  * Cao et al. [2013] W. Cao, J. Sun, and Z. Xu.  Fast image deconvolution using closed-form thresholding formulas of ${L}_{q}$ ($q=\tfrac{1}{2}$, $\tfrac{2}{3}$) regularization. _J. Vis. Commun. Image R._ , 24(1):31–41, 2013\.
  * Cartis et al. [2022] C. Cartis, N. I. M. Gould, and Ph. L. Toint.  _Evaluation Complexity of Algorithms for Nonconvex Optimization_.  Number 30 in MOS-SIAM Series on Optimization. SIAM, Philadelphia, USA, 2022.
  * Conn et al. [2000] 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, 2000\.
  * Dennis and Wolkowicz [1993] J. E. Dennis, Jr. and H. Wolkowicz.  Sizing and least-change secant methods.  _SIAM J. Numer. Anal._ , 30(5):1291–1314, 1993\.
  * Donoho [2006] D. Donoho.  Compressed sensing.  _IEEE T. Inform. Theory_ , 52(4):1289–1306, 2006\.
  * Duchi et al. [2011] J. Duchi, E. Hazan, and Y. Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  _J. Mach. Learn. Res._ , 12:2121–2159, 2011.
  * FitzHugh [1955] R. FitzHugh.  Mathematical models of threshold phenomena in the nerve membrane.  _B. Math. Biophys._ , 17(4):257–278, 1955.
  * Gilbert and Lemaréchal [1989] J.-C. Gilbert and C. Lemaréchal.  Some numerical experiments with variable-storage quasi-Newton algorithms.  _Math. Program._ , 45:407–435, 1989.
  * 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.
  * 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\.
  * Lotfi et al. [2020] S. Lotfi, T. Bonniot de Ruisselet, D. Orban, and A. Lodi.  Stochastic damped L-BFGS with controlled norm of the Hessian approximation.  2020\.  OPT2020 Conference on Optimization for Machine Learning.
  * 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\.
  * Nazareth [1995] J. L. Nazareth.  If quasi-Newton then why not quasi-Cauchy?  _SIAG/OPT Views-and-News_ , 6:11–14, 1995.
  * Orban et al. [2020] D. Orban, A. S. Siqueira, and contributors.  LinearOperators.jl.  https://github.com/JuliaSmoothOptimizers/LinearOperators.jl, September 2020.
  * Rockafellar and Wets [1998] R. Rockafellar and R. Wets.  _Variational Analysis_ , volume 317.  Springer Verlag, 1998.
  * Scheinberg and Tang [2016] K. Scheinberg and X. Tang.  Practical inexact proximal quasi-Newton method with global complexity analysis.  _Math. Program._ , (160):495–529, 2016.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _J. Roy. Statist. Soc. Ser. B_ , 58(1):267–288, 1996.
  * van der Pol [1926] B. van der Pol.  LXXXVIII. On “relaxation-oscillations”.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 2(11):978–992, 1926.
  * Zhu et al. [1999] M. Zhu, J. L. Nazareth, and H. Wolkowicz.  The quasi-Cauchy relation and diagonal updating.  _SIAM J. Optim._ , 9(4):1192–1204, 1999.