# An Inexact Projected Regularized Newton Method for Fused Zero-norms
Regularization Problems

Yuqia Wu111Department of Applied Mathematics, The Hong Kong Polytechnic
University, Hong Kong (yuqia.wu@connect.polyu.hk). Shaohua Pan222School of
Mathematics, South China University of Technology, Guangzhou
(shhpan@scut.edu.cn). Xiaoqi Yang333Department of Applied Mathematics, The
Hong Kong Polytechnic University, Hong Kong (mayangxq@polyu.edu.hk).

###### Abstract

We are concerned with structured $\ell_{0}$-norms regularization problems,
with a twice continuously differentiable loss function and a box constraint.
This class of problems have a wide range of applications in statistics,
machine learning and image processing. To the best of our knowledge, there is
no effective algorithm in the literature for solving them. In this paper, we
first obtain a polynomial-time algorithm to find a point in the proximal
mapping of the fused $\ell_{0}$-norms with a box constraint based on dynamic
programming principle. We then propose a hybrid algorithm of proximal gradient
method and inexact projected regularized Newton method to solve structured
$\ell_{0}$-norms regularization problems. The whole sequence generated by the
algorithm is shown to be convergent by virtue of a non-degeneracy condition, a
curvature condition and a Kurdyka-Łojasiewicz property. A superlinear
convergence rate of the iterates is established under a locally Hölderian
error bound condition on a second-order stationary point set, without
requiring the local optimality of the limit point. Finally, numerical
experiments are conducted to highlight the features of our considered model,
and the superiority of our proposed algorithm.

Keywords: fused $\ell_{0}$-norms regularization problems, polynomial-time
proximal gradient algorithm; inexact projected regularized Newton algorithm;
global convergence; superlinear convergence.

##  1 Introduction

Given a matrix $B\in\mathbb{R}^{p\times n}$, $\lambda_{1}>0$, $\lambda_{2}>0$,
$l\in\mathbb{R}^{n}_{-}$ and $u\in\mathbb{R}^{n}_{+}$, we are interested in
the following structured $\ell_{0}$-norms regularization problem with a box
constraint:

$\min_{x\in\mathbb{R}^{n}}\
F(x)\\!:=f(x)+\lambda_{1}\|Bx\|_{0}+\lambda_{2}\|x\|_{0}\ \ {\rm s.t.}\ \
l\leq x\leq u,$ (1)

where $f\\!:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a twice continuously
differentiable function, $\|\cdot\|_{0}$ denotes the $\ell_{0}$-norm (or
cardinality) function. This model encourages sparsity of both variable $x$ and
its linear transformation $Bx$. Throughout this paper, we define
$g(\cdot):=\lambda_{1}\|B\cdot\|_{0}+\lambda_{2}\|\cdot\|_{0}+\delta_{\Omega}(\cdot)$
and $\Omega:=\\{x\in\mathbb{R}^{n}\ \lvert\ l\leq x\leq u\\}$, where
$\delta_{\Omega}(\cdot)$ denotes the indicator function of $\Omega$.

###  1.1 Motivation

Given a data matrix $A\in\mathbb{R}^{m\times n}$ and its response
$b\in\mathbb{R}^{m}$, the usual regression model is to minimize
$f(x)\\!:=\\!h(Ax-b)$, where $h$ is a smooth function attaining its optimal
solution at the origin. For the case where
$h(\cdot)\\!=\\!\frac{1}{2}\|\cdot\|^{2}$, $f$ is the least-squares loss
function of the linear regression. It is known that one of the formulations
for finding a sparse vector while minimizing $f$ is the following $\ell_{0}$
regularization problem

$\min_{x\in\mathbb{R}^{n}}\ f(x)+\lambda_{2}\|x\|_{0},$ (2)

where the $\ell_{0}$-norm term shrinks some small coefficients to $0$ and
identifies a set of influential components. However, the $\ell_{0}$-norm
penalty only takes the sparsity of $x$ into consideration, but ignores its
spatial nature, which sometimes needs to be considered in real-world
applications. For example, in the context of image processing, the variables
often represent the pixels of images, which are correlated with their
neighboring ones. To recover the blurred images, Rudin et al. [35] took into
account the differences between adjacent variables and used the total
variation regularization, which penalizes the changes of the neighboring
pixels and hence encourages smoothness in the solution. Moreover, Land and
Friedman [21] studied the phoneme classification on TIMIT database (Acpistoc-
Phonetic Continuous Speech Corpus, NTIS, US Dept of Commerce), which consists
of 4509 32ms speech frames and each speech frame is represented by 512 samples
of 16 KHz rate. This database is collected from 437 male speakers. Every
speaker provided approximately two speech frames of each of five phonemes,
where the phonemes are “sh” as in “she”, “dcl” as in “dark”, “iy” as the vowel
in “she”, “aa” as the vowel in “dark”, and “ao” as the first vowel in “water”.
Since each phoneme is composed of a series of consecutively sampled points,
there is a high chance that each sampled point is close or identical to its
neighboring ones. For this reason, Land and Friedman [21] considered imposing
a fused penalty on the coefficients vector $x$, and proposed the following
problems with zero-order variable fusion and first-order variable fusion
respectively to train the classifier:

$\displaystyle\min_{x\in\mathbb{R}^{n}}\frac{1}{2}\|Ax-b\|^{2}+\lambda_{1}\|\widehat{B}x\|_{0},$
(3)
$\displaystyle\min_{x\in\mathbb{R}^{n}}\frac{1}{2}\|Ax-b\|^{2}+\lambda_{1}\|\widehat{B}x\|_{1},$
(4)

where $A\in\mathbb{R}^{m\times n}$ represents the phoneme data,
$b\in\mathbb{R}^{m}$ is the label vector,
$\widehat{B}\in\mathbb{R}^{(n\\!-\\!1)\times n}$ with $\widehat{B}_{ii}=1\
{\rm and}\ \widehat{B}_{i,i+1}=-1$ for all $i\in\\{1,\ldots,n\\!-\\!1\\}$ and
$\widehat{B}_{ij}=0$ otherwise. If
$f(\cdot)\\!=\\!\frac{1}{2}\|A\cdot-b\|^{2}$ and $B=\widehat{B}$, then we call
(1) a fused $\ell_{0}$-norms regularization problem with a box constraint.

Additionally taking the sparsity of $x$ into consideration, Tibshirani et al.
[41] proposed the fused Lasso, given by

$\min_{x\in\mathbb{R}^{n}}\
\frac{1}{2}\|Ax-b\|^{2}+\lambda_{1}\|\widehat{B}x\|_{1}+\lambda_{2}\|x\|_{1},$
(5)

and presented its nice statistical properties. Friedman et al. [15]
demonstrated that the proximal mapping of
$\lambda_{1}\|\widehat{B}x\|_{1}+\lambda_{2}\|x\|_{1}$ can be obtained through
a process, which is known as “prox-decomposition" later. Based on the
accessibility of this proximal mapping, various algorithms can efficiently
address model (5), see [24, 25, 23, 29]. In particular, Li et al. [23]
proposed a semismooth Newton augmented Lagrangian method (SSNAL) to solve the
dual of (5). The numerical results presented in their study indicate that
SSNAL is highly efficient.

It was claimed in [21] that both (3) and (4) perform well in signal
regression, but the zero-order fusion one produces simpler estimated
coefficient vectors. This observation suggests that model (1) with
$f=\frac{1}{2}\|A\cdot-b\|^{2}$ and $B=\widehat{B}$ may be able to effectively
find a simpler solution while performs well as the fused Lasso does. Compared
with regularization problems using $\ell_{0}$-norm, those using $\|Bx\|_{0}$
regularization remain less explored in terms of algorithm development.
According to [21], the global solution of (3) is unattainable. However, one of
its stationary points can be obtained. In fact, Jewell et al. [17] has
revealed by virtue of dynamic programming principle that a point in the
proximal mapping of $\lambda_{1}\|\widehat{B}\cdot\|_{0}$ can be exactly
determined within polynomial time, which allows one to use the well-known
proximal gradient (PG) method to find a stationary point of problem (3).
However, the highly nonconvex and nonsmooth nature of model (1) presents
significant challenges in computing the proximal mapping of $g$ when
$B=\widehat{B}$ and in developing effective optimization algorithms for
solving it. As far as we know, no specific algorithms have yet been designed
to solve these challenging problems.

Another motivation for this work comes from our previous research [45]. In
that work, we considered a model which replaces the $\ell_{0}$-norm in (2) by
$\|x\|_{q}^{q}$, where $q\in(0,1)$ and $\|x\|_{q}$ is the $\ell_{q}$ quasi-
norm of $x$, defined as $\|x\|_{q}=\big{(}\sum_{i=1}^{n}\lvert
x_{i}\rvert^{q}\big{)}^{1/q}$. For this class of nonconvex and nonsmooth
problems, we proposed a hybrid of PG and subspace regularized Newton methods
(HpgSRN), which in Newton step restricts the subproblem on a subspace within
which the objective function is smooth, and thus a regularized Newton method
can be applied. It is worth noting that the subspace is induced by the support
of the current iterate $x^{k}$. PG step is executed in every iteration, but it
does not necessarily run a Newton step unless a switch condition is satisfied.
A global convergence was established under a curvature condition and the
Kurdyka-Łojasiewicz (KL) property [3] of the objective function, and a
superlinear convergence rate was achieved under an additional local second-
order error bound condition. Due to the desirable convergence result and
numerical performance of HpgSRN, we aim to adapt a similar subspace
regularized Newton algorithm to solve (1), in which the subspace is induced by
the combined support of $Bx^{k}$ and $x^{k}$.

###  1.2 Related work

In recent years, many optimization algorithms have been well developed to
solve the $\ell_{0}$-norm regularization problems of the form (2), which
includes iterative hard thresholding [16, 10, 11, 27], the penalty
decomposition [28], the smoothing proximal gradient method [9], the
accelerated iterative hard thresholding [44] and NL0R [47]. Among all these
algorithms, NL0R is the only Newton-type method, which employs Newton method
to solve a series of stationary equations confined within the subspaces
identified by the support of the solution obtained by the proximal mapping of
$\lambda_{2}\|\cdot\|_{0}$.

The PG method is able to effectively process model (1) if the proximal mapping
of $g$ can be precisely computed. The PG method belongs to first-order
methods, which have a low computation cost and weak global convergence
conditions, but they achieve at most a linear convergence rate. On the other
hand, the Newton method has a faster convergence rate, but it can only be
applied to minimize sufficiently smooth objective functions. In recent years,
there have been active investigations into the Newton-type methods for
nonsmooth composite optimization problems of the form

$\min_{x\in\mathbb{R}^{n}}\ \Psi(x):=\psi(x)+\varphi(x),$ (6)

where
$\varphi:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}:=(-\infty,+\infty]$ is
proper lower semicontinuous, and $\psi(x)$ is twice continuously
differentiable on an open subset of $\mathbb{R}^{n}$ containing the domain of
$\varphi$. The proximal Newton-type method is able to address (6) with convex
$\varphi$ and convex or weakly convex $\psi$ (see [7, 22, 46, 30, 26]).
Another popular second-order method for solving (6) is to minimize the
forward-backward envelop (FBE) of $\Psi$, see [38, 40, 39, 1]. In particular,
for those $\Psi$ with the proximal mapping of $\varphi$ being available,
Themelis et al. [40] proposed an algorithm called ZeroFPR, based on the quasi-
Newton method, for minimizing the FBE of $\Psi$. They achieved the global
convergence of the iterate sequence by means of the KL property of the FBE and
its local superlinear rate under the Dennis-Moré condition and the strong
local minimality property of the limit point. An algorithm similar to ZeroFPR
but minimizing the Bregman FBE of $\Psi$ was proposed in [1], which achieves a
superlinear convergence rate without requiring the strong local minimality of
the limit point. For the case that $\psi$ is smooth and $\varphi$ admits a
computable proximal mapping, Bareilles et al. [5] proposed an algorithm,
alternating between a PG step and a Riemannian Newton method, which was proved
to have a quadratic convergence rate under a positive definiteness assumption
on the Riemannian Hessian at the limit point.

###  1.3 Main contributions

In this paper, we aim to design a hybrid of PG and inexact projected
regularized Newton methods (PGiPN) to solve the structured $\ell_{0}$-norms
regularization problem (1). Let $x^{k}\in\Omega$ be the current iterate. Our
method first runs a PG step with line search at $x^{k}$ to produce
$\overline{x}^{k}$ with

$\overline{x}^{k}\in{\rm
prox}_{\overline{\mu}_{k}^{-1}g}(x^{k}-\overline{\mu}_{k}^{-1}\nabla
f(x^{k})),$ (7)

where ${\rm prox}_{\overline{\mu}_{k}^{-1}g}(x)$ is the proximal mapping of
$g$ defined by (15), $\overline{\mu}_{k}>0$ is a constant such that $F$ gains
a sufficient decrease from $x^{k}$ to $\overline{x}^{k}$, and then judges
whether the iterate enters Newton step or not in terms of some switch
condition, which takes the following forms of structured stable supports:

${\rm supp}(x^{k})={\rm supp}(\overline{x}^{k})\ \ {\rm and}\ \ {\rm
supp}(Bx^{k})={\rm supp}(B\overline{x}^{k}).$ (8)

If this switch condition does not hold, we set $x^{k+1}=\overline{x}^{k}$ and
return to the PG step. Otherwise, due to the nature of $\ell_{0}$-norm, the
restriction of $\lambda_{1}\|Bx\|_{0}+\lambda_{2}\|x\|_{0}$ on the supports
${\rm supp}(Bx^{k})$ and ${\rm supp}(x^{k})$, i.e., $\lambda_{1}\|(Bx)_{{\rm
supp}(Bx^{k})}\|_{0}+\lambda_{2}\|x_{{\rm supp}(x^{k})}\|_{0}$, is a constant
near $x^{k}$ and does not provide any useful information at all and thus,
unlike dealing with the $\ell_{q}$ regularization problem in [45], we
introduce the following multifunction
$\Pi:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$:

$\displaystyle\Pi(z)$ $\displaystyle:=\left\\{x\in\Omega\ \lvert\ {\rm
supp}(x)\subset{\rm supp}(z),\ {\rm supp}(Bx)\subset{\rm supp}(Bz)\right\\}$
$\displaystyle=\left\\{x\in\Omega\ \lvert\ x_{[{\rm supp}(z)]^{c}}=0,\
(Bx)_{[{\rm supp}(Bz)]^{c}}=0\right\\},$ (9)

and consider the associated subproblem

$\min_{x\in\mathbb{R}^{n}}\,f(x)+\delta_{\Pi_{k}}(x)\ \ {\rm with}\ \
\Pi_{k}=\Pi(x^{k}).$ (10)

It is noted that the set $\Pi(x^{k})$ contains all the points whose supports
are a subset of the support of $x^{k}$ as well as the supports of their linear
transformation are a subset of the support of the linear transformation of
$x^{k}$. It is worth pointing out that the multifunction $\Pi$ is not closed
but closed-valued.

We will show that a stationary point of (10) is one for problem (1). Thus,
instead of a subspace regularized Newton step in [45], following the projected
Newton method in [7] and the proximal Newton method in [22, 46, 30] and [26],
our projected regularized Newton step minimizes the following second-order
approximation of (10) on $\Pi_{k}$:

$\mathop{\arg\min}_{x\in\mathbb{R}^{n}}\Big{\\{}\Theta_{k}(x)\\!:=f(x^{k})+\langle\nabla\\!f(x^{k}),x-\\!x^{k}\rangle+\frac{1}{2}\langle
x\\!-x^{k},G_{k}(x\\!-\\!x^{k})\rangle+\delta_{\Pi_{k}}(x)\Big{\\}},$ (11)

where $G_{k}$ is an approximation to the Hessian $\nabla^{2}\\!f(x^{k})$,
satisfying the following positive definiteness condition:

$G_{k}\succeq
b_{1}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}I,$ (12)

where $b_{1}>0$, $\sigma\in(0,\frac{1}{2})$ and $\overline{\mu}_{k}$ is the
one in (7). To cater for the practical computation, our Newton step seeks an
inexact solution $y^{k}$ of (11) satisfying

$\displaystyle\Theta_{k}(y)-\Theta_{k}(x^{k})\leq 0,$ (13) $\displaystyle{\rm
dist}(0,\partial\Theta_{k}(y))\leq\frac{\min\\{\overline{\mu}_{k}^{-1},1\\}}{2}\min\left\\{\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|,\|\overline{\mu}_{k}(x^{k}\\!-\overline{x}^{k})\|^{1+\varsigma}\right\\}$
(14)

with $\varsigma\in(\sigma,1]$. Setting the direction $d^{k}\\!:=y^{k}-x^{k}$,
a step size $\alpha_{k}\in(0,1]$ is found in the direction $d^{k}$ via
backtrackings, and set $x^{k+1}:=x^{k}+\alpha_{k}d^{k}$. To ensure the global
convergence, the next iterate still returns to the PG step. The details of the
algorithm are given in Section 3.

The main contributions of the paper are as follows:

$\bullet$ Based on dynamic programming principle, we develop a polynomial-time
($O(n^{3+o(1)})$) algorithm for seeking a point $\overline{x}_{k}$ in the
proximal mapping (7) of $g$ when $B=\widehat{B}$. This generalizes the
corresponding result in [17] for finding $\overline{x}^{k}$ in (7) from
$g(\cdot)=\lambda_{1}\|B\cdot\|_{0}$ to
$g(\cdot)=\lambda_{1}\|B\cdot\|_{0}+\lambda_{2}\|\cdot\|_{0}+\delta_{\Omega}(\cdot)$.
This also provides a PG algorithm for solving (1). We establish a uniform
lower bound on ${\rm prox}_{\mu^{-1}g}(x)$ for $x$ on a compact set and $\mu$
on a closed interval. This generalizes the corresponding results in [27] for
$\ell_{0}$-norm and in [45] for $\ell_{q}$-norm with $0\\!<\\!q\\!<\\!1$,
respectively.

$\bullet$ We design a hybrid algorithm (PGiPN) of PG and inexact projected
regularized Newton method to solve the structured $\ell_{0}$-norms
regularization problem (1), which includes the fused $\ell_{0}$-norms
regularization problem with a box constraint as a special case. We obtain the
global convergence of the algorithm by showing that the structured stable
supports (8) hold when the iteration number is sufficiently large. Moreover,
we establish a superlinear convergence rate under a Hölderian error bound on a
second-order stationary point set, without requiring the local optimality of
the limit point.

$\bullet$ The numerical experiments show that our PGiPN is more effective than
some existing algorithms in the literature in terms of solution quality and
efficiency.

The rest of the paper is organized as follows. In Section 2 we give some
preliminaries and characterizations of the optimality condition of model (1).
In Section 3, we give some results related to $g$, including a lower bound of
proximal mapping of $g$, and an algorithm for finding a point in the proximal
mapping of
$\lambda_{1}\|\widehat{B}x\|_{0}+\lambda_{2}\|x\|_{0}+\delta_{\Omega}(x)$. In
Section 4, we introduce our algorithm and show that it is well defined. In
Section 5 we present the convergence analysis of our algorithm. The
implementation scheme of our algorithm and the numerical experiments are
presented in Section 6.

##  2 Preliminaries

###  2.1 Notations

For any $x\in\mathbb{R}^{n}$ and $\epsilon>0$, $\mathbb{B}(x,\epsilon):=\\{z\
\lvert\ \|z-x\|\leq\epsilon\\}$ denotes the ball centered at $x$ with radius
$\epsilon$. Let ${\bf B}:=\mathbb{B}(0,1)$. For a closed and convex set
$\Xi\subset\mathbb{R}^{n}$, ${\rm proj}_{\Xi}(z)$ denotes the projection of a
point $z\in\mathbb{R}^{n}$ onto $\Xi$. We denote by $\mathcal{N}_{\Xi}(x)$ and
$\mathcal{T}_{\Xi}(x)$ the normal cone and tangent cone of $\Xi$ at $x$,
respectively. For a closed set $\Xi^{\prime}\subset\mathbb{R}^{n}$, ${\rm
dist}(z,\Xi^{\prime}):=\min_{x\in\Xi^{\prime}}\|x-z\|$. For $t\in\mathbb{R}$,
$t_{+}:=\max\\{t,0\\}$. Fix any two nonnegative integers $j<k$, define
$[j\\!:\\!k]:=\\{j,j\\!+\\!1,...,k\\}$ and $[k]:=[1\\!:\\!k]$. For an index
set $T\subset[n]$, write $T^{c}:=[n]\backslash T$ and $\lvert T\rvert$ is the
number of the elements of $T$. Given any $x\in\mathbb{R}^{n}$, ${\rm
supp}(x):=\\{i\in[n]\ \lvert\ x_{i}\neq 0\\},$ $\lvert
x\rvert_{\min}:=\min_{i\in{\rm supp}(x)}\lvert x_{i}\rvert$ and
$x_{T}\in\mathbb{R}^{\lvert T\rvert}$ is the vector consisting of those
$x_{j}$’s with $j\in T$, and $x_{j:k}:=x_{[j:k]}$. 1 and $I$ are the vector of
ones and the identity matrix, respectively, whose dimension is adaptive to the
context. Given a real symmetric matrix $H$, $\lambda_{\min}(H)$ denotes the
smallest eigenvalue of $H$, and $\|H\|_{2}$ is the spectral norm of $H$. For a
matrix $A\in\mathbb{R}^{m\times n}$ and $S\subset[m]$, $A_{S\cdot}$(resp.
$A_{\cdot T}$) denotes the matrix consisting of the rows (resp. the columns)
of $A$ whose indices correspond to $S$(resp. $T$). We write the range of $A$
by ${\rm Range}(A)=\\{Ax\ \lvert\ x\in\mathbb{R}^{n}\\}$ and the null space of
$A$ by ${\rm Null}(A)=\\{y\in\mathbb{R}^{n}\ \lvert\ Ay=0\\}$. For another
matrix $C\in\mathbb{R}^{m\times p}$, $[A\ C]\in\mathbb{R}^{m\times(n+p)}$ is
defined as a matrix composed of two matrices, $A$ and $C$, placed side by
side. For any $D\in\mathbb{R}^{p\times n}$, $[A;D]:=[A^{\top}\
D^{\top}]^{\top}$. For a proper lower semicontinuous function
$h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}:=\mathbb{R}\cup{\pm\infty}$,
its domain is denoted as ${\rm dom}h\\!:=\\{x\in\mathbb{R}^{n}\ \lvert\
h(x)<\infty\\}$, and the proximal mapping of $h$ is defined as

${\rm prox}_{\mu
h}(z):=\mathop{\arg\min}_{x\in\mathbb{R}^{n}}\frac{1}{2\mu}\|x-z\|^{2}+h(x),$
(15)

with $\mu>0$ and $z\in\mathbb{R}^{n}$.

###  2.2 Stationary conditions

Since problem (1) involves a box constraint and the structured
$\ell_{0}$-norms function is lower semicontinuous, the set of global optimal
solutions of model (1) is nonempty and compact. Moreover, by the continuity of
$\nabla^{2}\\!f$ and the compactness of $\Omega$, we have $\nabla\\!f$ is
Lipschitz continuous on $\Omega$, i.e., there exists $L_{1}>0$ such that

$\|\nabla f(x)-\nabla f(y)\|\leq L_{1}\|x-y\|\ \ {\rm for\ all}\
x,y\in\Omega.$ (16)

For a point $\overline{x}\in{\rm dom}h$, we denote the basic (or limiting)
subdifferential of $h$ at $\overline{x}$ by $\partial h(\overline{x})$ [34,
Definition 8.3]. Next, we introduce two classes of stationary points for the
general composite problem (6), which includes (1) as a special case.

###### Definition 2.1

A vector $x\in\mathbb{R}^{n}$ is called a stationary point of problem (6) if
$0\in\partial\Psi(x)$. A vector $x\in\mathbb{R}^{n}$ is called an
$L$-stationary point of problem (6) if there exists a constant $\mu>0$ such
that $x\in{\rm prox}_{\mu^{-1}\varphi}(x-\\!\mu^{-1}\nabla\psi(x))$.

Recall that $\Psi=\psi+\varphi$, where $\psi$ is twice continuously
differentiable and $\varphi$ is proper and lower semicontinuous. If in
addition $\varphi$ is assumed to be convex, then

$0\in\partial\Psi(x)\Leftrightarrow
0\in\mu(x-(x-\mu^{-1}\nabla\\!\psi(x)))+\partial\varphi(x)\Leftrightarrow
x={\rm prox}_{\mu^{-1}\varphi}(x\\!-\mu^{-1}\nabla\psi(x)),$

which means that for problem (6) the $L$-stationarity of a point $x$ is
equivalent to its stationarity. To extend this equivalence to a broader class
of functions, we recall the definition of prox-regularity, which acts as a
surrogate of convexity.

###### Definition 2.2

[34, Definition 13.27] A function
$h\\!:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ is prox-regular at a
point $\overline{x}\in{\rm dom}h$ for $\overline{v}\in\partial
h(\overline{x})$ if $h$ is locally lower semicontinuous at $\overline{x}$, and
there exist $r\geq 0$ and $\varepsilon>0$ such that $h(x^{\prime})\geq
h(x)+v^{\top}(x^{\prime}-x)-\frac{r}{2}\|x^{\prime}-x\|^{2}$ for all
$\|x^{\prime}-\overline{x}\|\leq\varepsilon$, whenever $v\in\partial h(x)$,
$\|v-\overline{v}\|<\varepsilon$, $\|x-\overline{x}\|<\varepsilon$ and
$h(x)<h(\overline{x})+\varepsilon$. If $h$ is prox-regular at $\overline{x}$
for all $\overline{v}\in\partial h(\overline{x})$, we say that $h$ is prox-
regular at $\overline{x}$.

The following proposition reveals that under the assumption of the prox-
regularity of $\varphi$, these two classes of stationary points for $\Psi$
coincide. The proof is similar to that in [45, Remark 2.5], and the details
are omitted here.

###### Proposition 2.1

If $\overline{x}$ is an $L$-stationary point of problem (6), then
$0\in\partial\Psi(\overline{x})$. If $\varphi$ is prox-regular at
$\overline{x}$ for $-\nabla\psi(\overline{x})$ and prox-bounded444For the
definition of prox-boundedness, see [34, Definitions 1.23]., the converse is
also true.

Next we focus on the stationary conditions of problem (1). Recall that the
multifunction $\Pi:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ defined in
(1.3) is closed-valued. The following lemma characterizes the subdifferential
of function $F$.

###### Lemma 2.1

Fix any $z\in\Omega$. The following statements are true.

(i)

$\partial F(z)=\nabla\\!f(z)+\partial
g(z)=\nabla\\!f(z)+\mathcal{N}_{\Pi(z)}(z)$.

(ii)

$0\in\nabla\\!f(x)+\mathcal{N}_{\Pi(z)}(x)$ implies that $0\in\partial F(x)$.

Proof: The first equality of part (i) follows by [34, Exercise 8.8], and the
second one uses [31, Lemma 2.2 (i)]. Next we consider part (ii). Let
$x\in\Pi(z)$. From the definition of $\Pi(\cdot)$, we have
$\Pi(x)\subset\Pi(z)$, which along with $x\in\Pi(x)$ implies that
$\mathcal{N}_{\Pi(z)}(x)\subset\mathcal{N}_{\Pi(x)}(x)$. Combining part (i),
we obtain the desired result. $\Box$

###### Remark 2.1

Lemma 2.1 (ii) provides a way to seek a stationary point of $F$. Indeed, for
any given $z\in\mathbb{R}^{n}$, if $x$ is a stationary point of problem
$\mathop{\arg\min}\\{f(y)\ \lvert\ y\in\Pi(z)\\}$, i.e., $0\in\nabla
f(x)+\mathcal{N}_{\Pi(z)}(x)$, then by Lemma 2.1 (ii) it necessarily satisfies
$0\in\partial F(x)$. This technique will be utilized in the design of our
algorithm. In particular, when obtaining a good estimate of the stationary
point, say $x^{k}$, we use a Newton step to minimize $f$ over the polyhedral
set $\Pi(x^{k})$, so as to enhance the speed of the algorithm.

###  2.3 Kurdyka-Łojasiewicz property

In this subsection, we introduce the concept of Kurdyka-Łojasiewicz (KL)
property. The KL property of a function plays an important role in the
convergence analysis of first-order algorithms for nonconvex and nonsmooth
optimization problems (see, e.g., [3, 4]). In this work, we will use it to
establish the global convergence property of our algorithm.

###### Definition 2.3

For any $\eta>0$, we denote by $\Upsilon_{\\!\eta}$ the set consisting of all
continuous concave $\varphi\\!:[0,\eta)\to\mathbb{R}_{+}$ that are
continuously differentiable on $(0,\eta)$ with $\varphi(0)=0$ and
$\varphi^{\prime}(s)>0$ for all $s\in(0,\eta)$. A proper function
$h\\!:\mathbb{R}^{n}\\!\to\overline{\mathbb{R}}$ is said to have the KL
property at $\overline{x}\in{\rm dom}\,\partial h$ if there exist
$\eta\in(0,\infty]$, a neighborhood $\mathcal{U}$ of $\overline{x}$ and a
function $\varphi\in\Upsilon_{\\!\eta}$ such that for all
$x\in\mathcal{U}\cap\big{[}h(\overline{x})<h<h(\overline{x})+\eta\big{]}$,
$\varphi^{\prime}(h(x)-h(\overline{x})){\rm dist}(0,\partial h(x))\geq 1.$ If
$h$ has the KL property at each point of ${\rm dom}\partial h$, then $h$ is
called a KL function.

##  3 Prox-regularity and proximal mapping of $g$

###  3.1 Prox-regularity of $g$

In this subsection, we aim at proving the prox-regularity of $g$, which
together with Proposition 2.1 and the prox-boundedness of $g$ indicates that
the set of stationary points of problem (1) coincides with that of its
$L$-stationary points.

We remark here that the prox-regularity of $g$ cannot be obtained from the
existing calculus of prox-regularity. In fact, it was revealed in [32, Theorem
3.2] that, for proper $f_{i},\ i=1,2$ with $f_{i}$ being prox-regular at
$\overline{x}$ for $v_{i}\in\partial f_{i}(\overline{x})$ and let
$v:=v_{1}+v_{2}$, and $f_{0}:=f_{1}+f_{2}$, a sufficient condition such that
$f_{0}$ is prox-regular at $\overline{x}$ for $v$ is

$w_{1}+w_{2}=0\ {\rm with}\
w_{i}\in\partial^{\infty}f_{i}(\overline{x})\Longrightarrow w_{i}=0,\ i=1,2,$
(17)

where $\partial^{\infty}$ denotes the horizon subdifferential [34, Definition
8.3]. We give a counter example to illustrate that the above constraint
qualification does not hold for $f_{i}:\mathbb{R}^{4}\rightarrow\mathbb{R}$
with $f_{1}=\|\widehat{B}\cdot\|_{0}$ and $f_{2}=\|\cdot\|_{0}$. Let
$\overline{x}=(0,0,0,1)^{\top}$. Then,

$\partial^{\infty}f_{1}(\overline{x})=\partial f_{1}(\overline{x})={\rm
Range}((\widehat{B}_{[2]\cdot})^{\top}),\
\partial^{\infty}f_{2}(\overline{x})=\partial f_{2}(\overline{x})={\rm
Range}((I_{[3]\cdot})^{\top}).$

By the expressions of $\partial^{\infty}f_{1}(\overline{x})$ and
$\partial^{\infty}f_{2}(\overline{x})$, it is immediate to check that the
constraint qualification in (17) does not hold. Next, we give our proof toward
the prox-regularity of $g$.

###### Lemma 3.1

The function $g$ is prox-regular on its domain $\Omega$. Consequently, the set
of stationary points of model (1) coincides with its set of $L$-stationary
points.

Proof: Fix any $\overline{x}\in\Omega$ and pick any $\overline{v}\in\partial
g(\overline{x})$. Let $\lambda:=\min\\{\lambda_{1},\lambda_{2}\\}$ and
$C:=[B;I]$. Pick any
$\varepsilon\in(0,\min\\{\lambda,\|\overline{v}\|,\frac{\lambda}{5\|\overline{v}\|}\\})$
such that for all $x\in\mathbb{B}(\overline{x},\varepsilon)$, ${\rm
supp}(Cx)\supset{\rm supp}(C\overline{x})$. Next we prove that

$g(x^{\prime})\geq g(x)+v^{\top}(x^{\prime}-x),\ {\rm for\ all\
}\|x^{\prime}-\overline{x}\|\leq\varepsilon,\ v\in\partial g(x),\
\|v-\overline{v}\|<\varepsilon\ {\rm and}\ x\in\Xi,$ (18)

where $\Xi:=\\{x\ |\ \|x-\overline{x}\|<\varepsilon,\
g(x)<g(\overline{x})+\varepsilon\\},$ which implies that $g$ is prox-regular
at $\overline{x}$ for $\overline{v}$.

We first claim that for each $x\in\Xi$, it holds that ${\rm supp}(Cx)={\rm
supp}(C\overline{x})$ and $x\in\Omega$. If fact, by the definition of
$\varepsilon$, ${\rm supp}(Cx)\supset{\rm supp}(C\overline{x})$. If ${\rm
supp}(Cx)\neq{\rm supp}(C\overline{x})$, we have $g(x)\geq
g(\overline{x})+\lambda>g(\overline{x})+\varepsilon$, which yields that
$x\notin\Xi$. Therefore, ${\rm supp}(Cx)={\rm supp}(C\overline{x})$. The fact
that $x\in\Xi$ implies $x\in\Omega$ is clear. Hence the claimed facts are
true.

Fix any $x\in\Xi$. Consider any
$x^{\prime}\in\mathbb{B}(\overline{x},\varepsilon)$. If
$x^{\prime}\notin\Omega$, since $g(x^{\prime})=\infty$, it is immediate to see
that (18) holds, so it suffices to consider
$x^{\prime}\in\mathbb{B}(\overline{x},\varepsilon)\cap\Omega$. Note that ${\rm
supp}(Cx^{\prime})\supset{\rm supp}(C\overline{x})={\rm supp}(Cx)$. If ${\rm
supp}(Cx^{\prime})\neq{\rm supp}(Cx)$, then $g(x^{\prime})\geq g(x)+\lambda$.
For any $v\in\partial g(x)$ with $v\in\mathbb{B}(\overline{v},\varepsilon)$,
$\|v\|\leq\|\overline{v}\|+\varepsilon\leq 2\|\overline{v}\|$, which along
with $\|x^{\prime}-x\|\leq\|x^{\prime}-\overline{x}\|+\|x-\overline{x}\|\leq
2\varepsilon$ implies that

$g(x^{\prime})-g(x)-v^{\top}(x^{\prime}-x)\geq\lambda-\|v\|\|x^{\prime}-x\|\geq\lambda-4\|\overline{v}\|\varepsilon>0.$

Equation (18) holds. Next we consider the case ${\rm supp}(Cx^{\prime})={\rm
supp}(Cx)$. Define

$\Pi^{1}(x):=\big{\\{}z\in\mathbb{R}^{n}\ \lvert\ (Bz)_{[{\rm
supp}(Bx)]^{c}}=0\big{\\}},\ \Pi^{2}(x):=\big{\\{}z\in\mathbb{R}^{n}\ \lvert\
z_{[{\rm supp}(x)]^{c}}=0\big{\\}}.$

Clearly, $\Pi(x)=\Pi^{1}(x)\cap\Pi^{2}(x)\cap\Omega$ and
$\Pi_{1}(x),\Pi_{2}(x)$ and $\Omega$ are all polyhedral sets. By [33, Theorem
23.8], for any $v\in\mathcal{N}_{\Pi(x)}(x)=\partial g(x)$, there exist
$v_{1}\in\mathcal{N}_{\Pi^{1}(x)}(x)$, $v_{2}\in\mathcal{N}_{\Pi^{2}(x)}(x)$
and $v_{3}\in\mathcal{N}_{\Omega}(x)$ such that $v=v_{1}+v_{2}+v_{3}$. Then,

$\displaystyle g(x^{\prime})-g(x)-$ $\displaystyle
v^{\top}(x^{\prime}-x)=\lambda_{1}\|Bx^{\prime}\|_{0}-\lambda_{1}\|Bx\|_{0}-v_{1}^{\top}(x^{\prime}-x)$
$\displaystyle+\lambda_{2}\|x^{\prime}\|_{0}-\lambda_{2}\|x\|_{0}-v_{2}^{\top}(x^{\prime}-x)-v_{3}^{\top}(x^{\prime}-x)\geq
0,$

where the inequality follows from
$\lambda_{1}\|Bx^{\prime}\|_{0}-\lambda_{1}\|Bx\|_{0}=0,v_{1}^{\top}(x^{\prime}-x)=0$,
$\lambda_{2}\|x^{\prime}\|_{0}-\lambda_{2}\|x\|_{0}=0,v_{2}^{\top}(x^{\prime}-x)=0$
and $v_{3}^{\top}(x^{\prime}-x)\leq 0$. Equation (18) is true. Thus, by the
arbitrariness of $\overline{x}\in\Omega$ and $\overline{v}\in\partial
g(\overline{x})$, we conclude that $g$ is prox-regular on set $\Omega$. $\Box$

###  3.2 Lower bound of the proximal mapping of $g$

Given $\lambda>0$ and $x\in\mathbb{R}^{n}$, for any $z\in{\rm
prox}_{\lambda\|\cdot\|_{0}}(x)$, it holds that if $|z_{i}|>0$, then
$|z_{i}|\geq\sqrt{2\lambda}$ [27, Lemma 3.3]. This indicates that $|z|_{\min}$
has a uniform lower bound on ${\rm prox}_{\lambda\|\cdot\|_{0}}(x)$. Such a
uniform lower bound is shown to hold for $\ell_{q}$-norm with
$0\\!<\\!q\\!<\\!1$ and played a crucial role in the convergence analysis of
the algorithms involving subspace Newton method (see [45]). Next, we show that
such a uniform lower bound exists for $g$.

###### Lemma 3.2

For any given compact set $\Xi\subset\mathbb{R}^{n}$ and constants
$0<\underline{\mu}<\overline{\mu}$, define

$\mathcal{Z}:={\textstyle\bigcup_{z\in\Xi,\mu\in[\underline{\mu},\overline{\mu}]}}\
{\rm prox}_{\mu^{-1}g}(z).$

Then, there exists $\nu>0$ (depending on $\Xi,\underline{\mu}$ and
$\overline{\mu}$) such that
$\inf_{u\in\mathcal{Z}\backslash\\{0\\}}\,\lvert[B;I]u\rvert_{\min}\geq\nu$.

Proof: Write $C:=[B;I]$. By invoking [6, Corollary 3] and the compactness of
$\Omega$, there exists $\kappa>0$ such that for all index set
$J\subset[n\\!+\\!p]$,

${\rm dist}(x,{\rm Null}(C_{J\cdot})\cap\Omega)\leq\kappa{\rm dist}(x,{\rm
Null}(C_{J\cdot}))\ \ {\rm for\ any\ }x\in\Omega.$ (19)

In addition, there exists $\sigma>0$ such that for any index set
$J\subset[n\\!+\\!p]$ with $\\{C_{j\cdot}\\}_{j\in J}$ being linearly
independent,

$\lambda_{\min}(C_{J\cdot}C_{J\cdot}^{\top})\geq\sigma.$ (20)

For any $z\in\Xi$ and $\mu\in[\underline{\mu},\overline{\mu}]$, define
$h_{z,\mu}(x):=\frac{\mu}{2}\|x-z\|^{2}$ for $x\in\mathbb{R}^{n}$. By the
compactness of $\Omega$, $[\underline{\mu},\overline{\mu}]$ and $\Xi$, there
exists $\delta_{0}\in(0,1)$ such that for all $z\in\Xi$,
$\mu\in[\underline{\mu},\overline{\mu}]$ and $x,y\in\Omega$ with
$\|x-y\|<\delta_{0}$,
$\overline{\mu}(\|x\|+\|y\|+2\|z\|)\|x-y\|<\lambda:=\min\\{\lambda_{1},\lambda_{2}\\}$,
and consequently,

$\displaystyle\lvert h_{z,\mu}(x)-h_{z,\mu}(y)\rvert$
$\displaystyle=\frac{\mu}{2}\lvert\langle
x-y,x+y-2z\rangle\rvert\leq\frac{\overline{\mu}}{2}(\|x\|+\|y\|+2\|z\|)\|x-y\|<\frac{\lambda}{2}.$
(21)

Now suppose that the conclusion does not hold. Then there is a sequence
$\\{\overline{z}^{k}\\}_{k\in\mathbb{N}}\subset\mathcal{Z}\backslash\\{0\\}$
such that $\lvert C\overline{z}^{k}\rvert_{\min}\leq\frac{1}{k}$ for all
$k\in\mathbb{N}$. Note that $C$ has a full column rank. We also have $\lvert
C\overline{z}^{k}\rvert_{\min}>0$ for each $k\in\mathbb{N}$. By the definition
of $\mathcal{Z}$, for each $k\in\mathbb{N}$, there exist $z^{k}\in\Xi$ and
$\mu_{k}\in[\underline{\mu},\overline{\mu}]$ such that
$\overline{z}^{k}\in{\rm prox}_{\mu_{k}^{-1}g}(z^{k})$. Since $\lvert
C\overline{z}^{k}\rvert_{\min}\in\big{(}0,\frac{1}{k}\big{)}$ for all
$k\in\mathbb{N}$, there exist $\mathcal{K}\subset\mathbb{N}$ and an index
$i\in[n\\!+\\!p]$ such that

$0<\lvert(C\overline{z}^{k})_{i}\rvert=\lvert
C\overline{z}^{k}\rvert_{\min}<\frac{\delta_{0}\sigma}{\kappa\|C\|_{2}}\quad{\rm
for\ each}\ k\in\mathcal{K},$ (22)

where $\kappa$ and $\sigma$ are the ones appearing in (19) and (20),
respectively. Fix any $k\in\mathcal{K}$. Write
$Q_{k}:=[n\\!+\\!p]\backslash{\rm supp}(C\overline{z}^{k})$ and choose
$J_{k}\subset Q_{k}$ such that the rows of $C_{J_{k}\cdot}$ form a basis of
those of $C_{Q_{k}\cdot}$. Let $\widehat{J}_{k}:=J_{k}\cup\\{i\\}$. If
$J_{k}=\emptyset$, then $C_{\widehat{J}_{k}\cdot}$ has a full row rank. If
$J_{k}\neq\emptyset$, then $C_{J_{k}\cdot}\overline{z}^{k}=0$, which implies
that $C_{\widehat{J}_{k}\cdot}$ also has a full row rank (if not, $C_{i\cdot}$
is a linear combination of the rows of $C_{J_{k}\cdot}$, which along with
$C_{J_{k}\cdot}\overline{z}^{k}=0$ implies that
$C_{i\cdot}\overline{z}^{k}=0$, contradicting to
$\lvert(C\overline{z}^{k})_{i}\rvert=\lvert C\overline{z}^{k}\rvert_{\rm
min}>0$). Let $\widetilde{z}^{k}:={\rm proj}_{{\rm
Null}(C_{\widehat{J}_{k}\cdot})}(\overline{z}^{k})$. Then,
$C_{\widehat{J}_{k}\cdot}\widetilde{z}^{k}=0$ and
$(\overline{z}^{k}-\widetilde{z}^{k})\in{\rm
Range}\,(C_{\widehat{J}_{k}\cdot}^{\top})$. The latter means that there exists
$\xi^{k}\in\mathbb{R}^{\lvert\widehat{J}_{k}\rvert}$ such that
$\overline{z}^{k}-\widetilde{z}^{k}=C_{\widehat{J}_{k}\cdot}^{\top}\xi^{k}$.
Since $C_{\widehat{J}_{k}\cdot}$ has a full row rank and
$\|C_{\widehat{J}_{k}\cdot}\overline{z}^{k}\|=\lvert(C\overline{z}^{k})_{i}\rvert$,
we have

$\lvert(C\overline{z}^{k})_{i}\rvert=\|C_{\widehat{J}_{k}\cdot}\overline{z}^{k}-C_{\widehat{J}_{k}\cdot}\widetilde{z}^{k}\|=\|C_{\widehat{J}_{k}\cdot}C_{\widehat{J}_{k}\cdot}^{\top}\xi^{k}\|\geq\sigma\|\xi^{k}\|,$
(23)

where the last inequality is due to (20). Combining (23) with (22) yields
$\|\xi^{k}\|<\kappa^{-1}\|C\|_{2}^{-1}\delta_{0}$. Therefore,

$\|\overline{z}^{k}-\widetilde{z}^{k}\|=\|C_{\widehat{J}_{k}\cdot}^{\top}\xi^{k}\|\leq\|C_{\widehat{J}_{k}\cdot}\|_{2}\|\xi^{k}\|\leq\|C\|_{2}\|\xi^{k}\|<\kappa^{-1}\delta_{0}.$
(24)

Let $\widehat{z}^{k}:={\rm proj}_{{\rm
Null}(C_{\widehat{J}_{k}\cdot})\cap\Omega}(\overline{z}^{k})$. From (19) and
(24), it follows that

$\|\overline{z}^{k}-\widehat{z}^{k}\|={\rm dist}(\overline{z}^{k},{\rm
Null}(C_{\widehat{J}_{k}\cdot})\cap\Omega)\leq\kappa{\rm
dist}(\overline{z}^{k},{\rm
Null}(C_{\widehat{J}_{k}\cdot}))=\kappa\|\overline{z}^{k}-\widetilde{z}^{k}\|<\delta_{0}.$
(25)

Note that $\widehat{z}^{k},\overline{z}^{k}\in\Omega$. From (25) and (21), it
follows that

$\lvert
h_{z^{k},\mu_{k}}(\widehat{z}^{k})-h_{z^{k},\mu_{k}}(\overline{z}^{k})\rvert<\frac{\lambda}{2}.$
(26)

Next we claim that ${\rm supp}(C\widehat{z}^{k})\cup\\{i\\}\subset{\rm
supp}(C\overline{z}^{k})$. Indeed, since the rows of
$C_{\widehat{J}_{k}\cdot}$ form a basis of those of
$C_{[Q_{k}\cup\\{i\\}]\cdot}$ and $C_{\widehat{J}_{k}\cdot}\widehat{z}^{k}=0$,
$C_{[Q_{k}\cup\\{i\\}]\cdot}\widehat{z}^{k}=0$. Then, ${\rm
supp}(C_{[Q_{k}\cup\\{i\\}]\cdot}\widehat{z}^{k})\cup\\{i\\}={\rm
supp}(C_{[Q_{k}\cup\\{i\\}]\cdot}\overline{z}^{k})$. Since all the entries of
$C_{[Q_{k}\cup\\{i\\}]^{c}\cdot}\overline{z}^{k}$ are nonzero, it holds that
${\rm supp}(C_{[Q_{k}\cup\\{i\\}]^{c}\cdot}\widehat{z}^{k})\subset{\rm
supp}(C_{[Q_{k}\cup\\{i\\}]^{c}\cdot}\overline{z}^{k})$, which implies that
${\rm supp}(C\widehat{z}^{k})\cup\\{i\\}\subset{\rm supp}(C\overline{z}^{k})$.
Thus, the claimed inclusion follows, which implies that
$g(\overline{z}^{k})-g(\widehat{z}^{k})\geq\lambda$. This together with (26)
yields
$h_{z^{k},\mu_{k}}(\overline{z}^{k})+g(\overline{z}^{k})-(h_{z^{k},\mu_{k}}(\widehat{z}^{k})+g(\widehat{z}^{k}))\geq\lambda-\frac{\lambda}{2}=\frac{\lambda}{2}$,
contradicting to $\overline{z}^{k}\in{\rm prox}_{\mu_{k}^{-1}g}(z^{k})$. The
proof is completed. $\Box$

The result of Lemma 3.2 will be utilized in Proposition 4.1 to justify the
fact that the sequences $\\{\lvert
B\overline{x}^{k}\rvert_{\min}\\}_{k\in\mathbb{N}}$ and
$\\{\lvert\overline{x}^{k}\rvert_{\min}\\}_{k\in\mathbb{N}}$ are uniformly
lower bounded, where $\overline{x}^{k}$ is obtained in (7) (or (37) below).
This is a crucial aspect in proving the stability of ${\rm supp}(x^{k})$ and
${\rm supp}(Bx^{k})$ when $k$ is sufficiently large.

###  3.3 Proximal mapping of a fused $\ell_{0}$-norms function with a box
constraint

Using the idea of [20], Jewell et al. [17] presented a polynomial-time
algorithm for computing the proximal mapping of the fused $\ell_{0}$-norm
$\lambda_{1}\|\widehat{B}\cdot\|_{0}$, where
$\widehat{B}x=(x_{1}-x_{2};...;x_{n-1}-x_{n})$ for any $x\in\mathbb{R}^{n}$.
We extend the result of [17] for computing the proximal mapping of the fused
$\ell_{0}$-norms
$\lambda_{1}\|\widehat{B}\cdot\|_{0}+\lambda_{2}\|\cdot\|_{0}+\delta_{\Omega}(\cdot)$,
i.e., for any given $z\in\mathbb{R}^{n}$, seeking a global optimal solution of
the problem

$\min_{x\in\mathbb{R}^{n}}\,h(x;z):=\frac{1}{2}\|x-z\|^{2}+\lambda_{1}\|\widehat{B}x\|_{0}+\lambda_{2}\|x\|_{0}+\delta_{\Omega}(x).$
(27)

To simplify the deduction, for each $i\in[n]$, we define
$\omega_{i}:\mathbb{R}\rightarrow\overline{\mathbb{R}}$ by
$\omega_{i}(\alpha):=\lambda_{2}\lvert\alpha\rvert_{0}+\delta_{[l_{i},u_{i}]}(\alpha)$.
It is clear that for all $x\in\mathbb{R}^{n}$,
$\lambda_{2}\|x\|_{0}+\delta_{\Omega}(x)=\sum_{i=1}^{n}\omega_{i}(x_{i})$. Let
$H(0):=-\lambda_{1}$, and for each $s\in[n]$, define

$H(s):=\min_{y\in\mathbb{R}^{s}}h_{s}(y;z_{1:s}),\ {\rm where}\
h_{s}(y;z_{1:s}):=\frac{1}{2}\|y-z_{1:s}\|^{2}+\lambda_{1}\|\widehat{B}_{\cdot[s]}y\|_{0}+\sum_{j=1}^{s}\omega_{j}(y_{j}).$
(28)

It is immediate to see that $H(n)$ is the optimal value to (27). For each
$s\in[n]$, define function
$P_{s}\\!:[0\\!:\\!s\\!-\\!1]\times\mathbb{R}\to\overline{\mathbb{R}}$ by

$\displaystyle P_{s}(i,\alpha):=H(i)+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s}\|^{2}+\sum_{j=i+1}^{s}\omega_{j}(\alpha)+\lambda_{1}.$ (29)

For any given $y\in\mathbb{R}^{s}$, if $i$ is the largest integer in
$[0\\!:\\!s\\!-\\!1]$ such that $y_{i}\neq y_{i+1}$, and
$y_{i+1}=y_{i+2}=...=y_{s}=\alpha$, then $y=(y_{1:i};\alpha{\bf 1})$ and

$\displaystyle h_{s}(y;z_{1:s})$
$\displaystyle\\!=\\!\frac{1}{2}\|y_{1:i}\\!-\\!z_{1:i}\|^{2}\\!+\\!\lambda_{1}\|\widehat{B}_{\cdot[i]}y_{1:i}\|_{0}\\!+\\!\sum_{j=1}^{i}\omega_{j}(y_{j})\\!+\\!\frac{1}{2}\|y_{i+1:s}\\!-\\!z_{i+1:s}\|^{2}\\!+\\!\sum_{j=i+1}^{s}\omega_{j}(y_{j})+\lambda_{1}$
$\displaystyle=h_{i}(y_{1:i};z_{1:i})+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s}\|^{2}+\sum_{j=i+1}^{s}\omega_{j}(\alpha)+\lambda_{1}.$

If $y_{1:i}$ is optimal to
$\min_{y^{\prime}\in\mathbb{R}^{i}}h_{i}(y^{\prime};z_{1:i})$, then
$H(i)=h_{i}(y_{1:i};z_{1:i})$, which by the definitions of $P_{s}$ and $h_{s}$
yields that $P_{s}(i,\alpha)=h_{s}(y;z_{1:s}).$ In the following lemma, we
prove that the optimal value of
$\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)$ is equal to $H(s)$, by
which we give characterization to an optimal solution of
$h_{s}(\cdot;z_{1:s})$.

###### Lemma 3.3

Fix any $s\in[n]$. The following statements are true.

(i)

$H(s)=\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)$.

(ii)

Assume that
$(i_{s}^{*},\alpha_{s}^{*})\in\mathop{\arg\min}_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)$.
Then $y^{*}=(y^{*}_{1:i_{s}^{*}};\alpha_{s}^{*}{\bf 1})$ is a global solution
of $\min_{y\in\mathbb{R}^{s}}h_{s}(y;z_{1:s})$ with
$y_{1:i_{s}^{*}}^{*}\in\mathop{\arg\min}_{v\in\mathbb{R}^{i_{s}^{*}}}h_{i_{s}^{*}}(v;z_{1:i_{s}^{*}})$.

Proof: (i) Let $y^{*}$ be an optimal solution to problem (28). If
$y^{*}_{i}=y^{*}_{j}$ for all $i,j\in[s]$, let $i_{s}^{*}=0$; otherwise, let
$i^{*}_{s}$ be the largest integer such that $y^{*}_{i^{*}_{s}}\neq
y^{*}_{i^{*}_{s}+1}$. Set $\alpha_{s}^{*}=y^{*}_{i^{*}_{s}+1}$. If
$i^{*}_{s}\neq 0$, from the definition of $H(\cdot)$,
$h_{i^{*}_{s}}(y^{*}_{1:i^{*}_{s}};z_{1:i^{*}_{s}})\geq H(i^{*}_{s})$, which
implies that

$\displaystyle\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)\leq
H(i^{*}_{s})+\frac{1}{2}\|\alpha^{*}_{s}{\bf
1}-z_{i_{s}^{*}+1:s}\|^{2}+\sum_{j=i^{*}_{s}+1}^{s}\omega_{j}(\alpha^{*}_{s})+\lambda_{1}$
$\displaystyle\leq
h_{i_{s}^{*}}(y^{*}_{1:i^{*}_{s}};z_{1:i^{*}_{s}})+\frac{1}{2}\|y^{*}_{i^{*}_{s}+1:s}-z_{i_{s}^{*}+1:s}\|^{2}+\sum_{j=i^{*}_{s}+1}^{s}\omega_{j}(y^{*}_{j})+\lambda_{1}=h_{s}(y^{*};z_{1:s})=H(s),$

where the first equality holds by $y^{*}_{i^{*}_{s}+1}\neq y^{*}_{i^{*}_{s}}$
and the expression of $h_{s}(y^{*};z_{1:s})$. If $i^{*}_{s}=0$,

$\displaystyle\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)\leq
H(0)+\frac{1}{2}\|y^{*}-z_{1:s}\|^{2}+\sum_{j=1}^{s}\omega_{j}(y^{*}_{j})+\lambda_{1}=H(s).$

Therefore, $\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)\leq H(s)$
holds. On the other hand, let $(i_{s}^{*},\alpha^{*}_{s})$ be an optimal
solution to $\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)$. If
$i_{s}^{*}\neq 0$, let $y^{*}\in\mathbb{R}^{s}$ be such that
$y^{*}_{1:i_{s}^{*}}\in\mathop{\arg\min}_{v\in\mathbb{R}^{i_{s}^{*}}}h_{i_{s}^{*}}(v;z_{1:i_{s}^{*}})$
and $y^{*}_{i_{s}^{*}+1:s}=\alpha_{s}^{*}{\bf 1}$. Then, it is clear that

$\displaystyle H(s)$ $\displaystyle\leq h_{s}(y^{*};z_{1:s})\leq
h_{i_{s}^{*}}(y^{*}_{1:i_{s}^{*}};z_{1:i_{s}^{*}})+\frac{1}{2}\|y^{*}_{i_{s}^{*}+1:s}-z_{i_{s}^{*}+1:s}\|^{2}+\sum_{j=i_{s}^{*}+1}^{s}\omega_{j}(y^{*}_{j})+\lambda_{1}$
$\displaystyle=H(i_{s}^{*})+\frac{1}{2}\|\alpha_{s}^{*}{\bf
1}-z_{i_{s}^{*}+1:s}\|^{2}+\sum_{j=i_{s}^{*}+1}^{s}\omega_{j}(\alpha_{s}^{*})+\lambda_{1}=\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha).$

If $i_{s}^{*}=0$, let $y^{*}=\alpha_{s}^{*}{\bf 1}$. We have

$H(s)\leq
h_{s}(y^{*};z_{1:s})=H(0)+\frac{1}{2}\|y^{*}-z_{1:s}\|^{2}+\sum_{j=1}^{s}\omega_{j}(\alpha_{s}^{*})+\lambda_{1}=\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha).$

Therefore, $H(s)\leq\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)$.
These two inequalities imply the result.

(ii) If $i_{s}^{*}\neq 0$, by part (i) and the definitions of $\alpha_{s}^{*}$
and $i_{s}^{*}$,

$\displaystyle H(s)$
$\displaystyle=\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)=H(i_{s}^{*})+\frac{1}{2}\|\alpha_{s}^{*}{\bf
1}-z_{i_{s}^{*}+1:s}\|^{2}+\sum_{j=i_{s}^{*}+1}^{s}\omega_{j}(\alpha_{s}^{*})+\lambda_{1}$
$\displaystyle=h_{i^{*}_{s}}(y^{*}_{1:i_{s}^{*}};z_{1:i_{s}^{*}})+\frac{1}{2}\|y^{*}_{i_{s}^{*}+1:s}-z_{i_{s}^{*}+1:s}\|^{2}+\sum_{j=i_{s}^{*}+1}^{s}\omega_{j}(y^{*}_{j})+\lambda_{1}\geq
h_{s}(y^{*},z_{1:s}),$

where the last inequality follows by the definition of $h_{s}(\cdot,z_{1:s})$.
If $i_{s}^{*}=0$,

$H(s)=\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha)=H(0)+\frac{1}{2}\|y^{*}-z_{1:s}\|^{2}+\sum_{j=1}^{s}\omega_{j}(y^{*}_{j})+\lambda_{1}=h_{s}(y^{*};z_{1:s}),$

Therefore, $H(s)\geq h_{s}(y^{*};z_{1:s}).$ Along with the definition of
$H(s)$, $H(s)=h_{s}(y^{*};z_{1:s})$. $\Box$

From Lemma 3.3 (i), the nonconvex nonsmooth problem (27) can be recast as a
mixed-integer programming with objective function given in (29). Lemma 3.3
(ii) suggests a recursive method to obtain an optimal solution to (27). In
fact, by setting $s=n$, we can obtain that there exists an optimal solution to
(27), says $x^{*}$, such that $x^{*}_{i_{n}^{*}+1:n}=\alpha_{n}^{*}{\bf 1}$,
and
$x^{*}_{1:i_{n}^{*}}\in\mathop{\arg\min}_{v\in\mathbb{R}^{i_{n}^{*}}}h_{i_{n}^{*}}(v;z_{1:i_{n}^{*}})$.
Next, by setting $s=i_{n}^{*}$, we are able to obtain the expression of
$x^{*}_{i_{i_{n}^{*}}^{*}+1:i_{n}^{*}}$. Repeating this loop backward until
$s=0$, we can obtain the full expression of an optimal solution to (27). The
outline of computing ${\rm
prox}_{\lambda_{1}\|\widehat{B}\cdot\|_{0}+\omega(\cdot)}(z)$ is shown as
follows.

$\displaystyle{\rm Set\ the\ current\ changepoint}\ s=n.$
$\displaystyle{\bf\rm While}\ s>0\ {\bf\rm do}$ $\displaystyle\quad\quad{\rm
Find}\
(i_{s}^{*},\alpha_{s}^{*})\in\mathop{\arg\min}_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha).$
(30) $\displaystyle\quad\quad{\rm Let}\
x^{*}_{i_{s}^{*}+1:s}=\alpha_{s}^{*}{\bf 1}\ {\rm and}\ s\leftarrow
i_{s}^{*}.$ $\displaystyle{\rm End}$

To obtain an optimal solution to (27), the remaining issue is how to execute
the first line in while loop of (3.3), or in other words, for any given
$s\in[n]$, how to find
$(i^{*}_{s},\alpha_{s}^{*})\in\mathbb{N}\times\mathbb{R}$ appearing in Lemma
3.3 (ii). The following proposition provides some preparations.

###### Proposition 3.1

For each $s\in[n]$, let
$P_{s}^{*}(\alpha):=\min_{i\in[0:s-1]}P_{s}(i,\alpha)$.

(i)

For all $\alpha\in\mathbb{R}$,

$P_{s}^{*}(\alpha)=\left\\{\begin{array}[]{ll}\frac{1}{2}(\alpha-
z_{1})^{2}+\omega_{1}(\alpha)&{\rm if}\ s=1,\\\
\min\Big{\\{}P_{s-1}^{*}(\alpha),\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})\\!+\\!\lambda_{1}\Big{\\}}\\!+\\!\frac{1}{2}(\alpha\\!-\\!z_{s})^{2}\\!+\\!\omega_{s}(\alpha)&{\rm
if}\ s\in[2\\!:\\!n].\end{array}\right.$

(ii)

Let $\mathcal{R}_{1}^{0}\\!:=\mathbb{R}$, and
$\mathcal{R}_{s}^{i}\\!:=\mathcal{R}_{s-1}^{i}\cap(\mathcal{R}_{s}^{s-1})^{c}$
for all $s\in[2\\!:\\!n]$ and $i\in[0\\!:\\!s\\!-\\!2]$, where

$\mathcal{R}_{s}^{s-1}:=\left\\{\alpha\in\mathbb{R}\ \lvert\
P_{s-1}^{*}(\alpha)\geq\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})+\lambda_{1}\right\\}.$
(31)

(a) For each $s\in[2\\!:\\!n]$,
$\bigcup_{i\in[0:s-1]}\mathcal{R}_{s}^{i}=\mathbb{R}$ and
$\mathcal{R}_{s}^{i}\cap\mathcal{R}_{s}^{j}=\emptyset$ for any $i\neq
j\in[0\\!:\\!s\\!-\\!1]$.
(b) For each $s\in[n]$ and $i\in[0\\!:\\!s\\!-\\!1]$,
$P_{s}^{*}(\alpha)=P_{s}(i,\alpha)$ when $\alpha\in\mathcal{R}_{s}^{i}$.

Proof: (i) Note that
$P_{1}^{*}(\alpha)=P_{1}(0,\alpha)=H(0)+\frac{1}{2}(\alpha-
z_{1})^{2}+\omega_{1}(\alpha)+\lambda_{1}=\frac{1}{2}(\alpha-
z_{1})^{2}+\omega_{1}(\alpha).$ Now fix any $s\in[2:n]$. By the definition of
$P_{s}^{*}$, for any $\alpha\in\mathbb{R}$,

$P_{s}^{*}(\alpha)=\min_{i\in[0:s-1]}P_{s}(i,\alpha)=\min\Big{\\{}\min_{i\in[0:s-2]}P_{s}(i,\alpha),\,P_{s}(s\\!-\\!1,\alpha)\Big{\\}}.$
(32)

From the definition of $P_{s}$ in (29), for each $i\in[0\\!:\\!s\\!-\\!2]$ and
$\alpha\in\mathbb{R}$, it holds that

$\displaystyle P_{s}(i,\alpha)$ $\displaystyle=H(i)+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s}\|^{2}+\sum_{j=i+1}^{s}\omega_{j}(\alpha)+\lambda_{1}$
$\displaystyle=H(i)+\frac{1}{2}\|\alpha{\bf
1}\\!-\\!z_{i+1:s-1}\|^{2}+\sum_{j=i+1}^{s-1}\omega_{j}(\alpha)+\lambda_{1}+\frac{1}{2}(\alpha\\!-\\!z_{s})^{2}+\omega_{s}(\alpha)$
$\displaystyle=P_{s-1}(i,\alpha)+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha),$

while for any $\alpha\in\mathbb{R}$,
$P_{s}(s\\!-\\!1,\alpha)=H(s\\!-\\!1)+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha)+\lambda_{1}.$ By combining the last two
equalities with (32), we immediately obtain that

$\displaystyle P_{s}^{*}(\alpha)$
$\displaystyle=\min\Big{\\{}\min_{i\in[0:s-2]}P_{s-1}(i,\alpha),\,H(s\\!-\\!1)+\lambda_{1}\Big{\\}}+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha)$ (33)
$\displaystyle=\min\Big{\\{}P_{s-1}^{*}(\alpha),\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})+\lambda_{1}\Big{\\}}+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha),$

where the last equality follows by Lemma 3.3 (i). Thus, we get the desired
result.

(ii) We first prove (a) by induction. When $s=2$, since
$\mathcal{R}_{1}^{0}=\mathbb{R}$ and
$\mathcal{R}_{2}^{0}=\mathcal{R}_{1}^{0}\cap(\mathcal{R}_{2}^{1})^{c}$, we
have $\mathcal{R}_{2}^{0}\cup\mathcal{R}_{2}^{1}=\mathbb{R}$ and
$\mathcal{R}_{2}^{0}\cap\mathcal{R}_{2}^{1}=\emptyset.$ Assume that the result
holds when $s=j$ for some $j\in[2\\!:\\!n\\!-\\!1]$. We consider the case
$s=j\\!+\\!1$. Since
$\mathcal{R}_{s}^{i}\\!:=\mathcal{R}_{s-1}^{i}\cap(\mathcal{R}_{s}^{s-1})^{c}$
for all $i\in[0\\!:\\!s\\!-\\!2]$ and
$\bigcup_{i\in[0:s-2]}\mathcal{R}_{s-1}^{i}=\mathbb{R}$, it holds that

$\displaystyle{\textstyle\bigcup_{i\in[0:s-1]}}\mathcal{R}_{s}^{i}$
$\displaystyle=\big{[}{\textstyle\bigcup_{i\in[0:s-2]}(\mathcal{R}_{s-1}^{i}\cap(\mathcal{R}_{s}^{s-1})^{c})}\big{]}\cup\mathcal{R}_{s}^{s-1}=\big{(}\mathbb{R}\cap(\mathcal{R}_{s}^{s-1})^{c}\big{)}\cup\mathcal{R}_{s}^{s-1}=\mathbb{R}.$

Thus we obtain the first part of (a) by deduction. For any
$i\in[0\\!:\\!s\\!-\\!2]$, by definition,
$\mathcal{R}_{s}^{i}\cap\mathcal{R}_{s}^{s-1}=\emptyset$. It suffices to show
that $\mathcal{R}_{s}^{i}\cap\mathcal{R}_{s}^{j}=\emptyset$ for any $i\neq
j\in[0\\!:\\!s\\!-\\!2]$. By definition,

$\mathcal{R}_{s}^{i}\cap\mathcal{R}_{s}^{j}=\big{[}\mathcal{R}_{s-1}^{i}\cap(\mathcal{R}^{s-1}_{s})^{c}\big{]}\cap\big{[}\mathcal{R}_{s-1}^{j}\cap(\mathcal{R}^{s-1}_{s})^{c}\big{]}=\emptyset,$

where the last equality is using
$\mathcal{R}_{s-1}^{i}\cap\mathcal{R}_{s-1}^{j}=\emptyset.$ Thus, the second
part of (a) is obtained.

Next we prove (b). Since for any $\alpha\in\mathbb{R}=\mathcal{R}_{1}^{0}$,
$P_{1}^{*}(\alpha)=P_{1}(0,\alpha)$, the result holds for $s=1$. For
$s\in[2\\!:\\!n]$ and $i\\!=\\!s\\!-\\!1$, by the definition of
$\mathcal{R}_{s}^{s-1}$, for all $\alpha\in\mathcal{R}_{s}^{s-1}$,

$P_{s}^{*}(\alpha)=\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})+\lambda_{1}+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha)=P_{s}(s-1,\alpha),$

where the second equality is using Lemma 3.3 (i) and the definition of
$P_{s}$. Next we consider $s\in[2\\!:\\!n]$ and $i\in[0\\!:\\!s\\!-\\!2]$. We
argue by induction that $P_{s}^{*}(\alpha)=P_{s}(i,\alpha)$ when
$\alpha\in\mathcal{R}_{s}^{i}$. Indeed, when $s=2$, since
$\mathcal{R}_{2}^{0}=\mathcal{R}_{1}^{0}\cap(\mathcal{R}_{2}^{1})^{c}=(\mathcal{R}_{2}^{1})^{c}$,
for any $\alpha\in\mathcal{R}_{2}^{0}$, from (31) we have
$P_{1}^{*}(\alpha)<\min_{\alpha^{\prime}\in\mathbb{R}}P_{1}^{*}(\alpha^{\prime})+\lambda_{1}$,
which by part (i) implies that
$P_{2}^{*}(\alpha)=P_{1}^{*}(\alpha)+\frac{1}{2}(\alpha-
z_{2})^{2}+\omega_{2}(\alpha)=P_{1}(0,\alpha)+\frac{1}{2}(\alpha-
z_{2})^{2}+\omega_{2}(\alpha)=P_{2}(0,\alpha)$. Assume that the result holds
when $s=j$ for some $j\in[2\\!:\\!n\\!-\\!1]$. We consider the case for
$s=j\\!+\\!1$. For any $i\in[0\\!:\\!s\\!-\\!2]$, by definition,
$\mathcal{R}_{s}^{i}=\mathcal{R}_{s-1}^{i}\cap(\mathcal{R}_{s}^{s-1})^{c}$.
Then, from (33) for any $\alpha\in\mathcal{R}_{s}^{i}$,

$\displaystyle P_{s}^{*}(\alpha)$
$\displaystyle=P_{s-1}^{*}(\alpha)+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha)=P_{s-1}(i,\alpha)+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha)$ $\displaystyle=H(i)+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s-1}\|^{2}+\sum_{j=i+1}^{s-1}w_{j}(\alpha)+\lambda_{1}+\frac{1}{2}(\alpha\\!-\\!z_{s})^{2}+\omega_{s}(\alpha)$
$\displaystyle=H(i)+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s}\|^{2}+\sum_{j=i+1}^{s}w_{j}(\alpha)+\lambda_{1}=P_{s}(i,\alpha),$

where the second equality is using $P_{s-1}^{*}(\alpha)=P_{s-1}(i,\alpha)$
implied by induction. Hence, the conclusion holds for $s=j+1$ and any
$i\in[0\\!:\\!s\\!-\\!2]$. The proof is completed. $\Box$

Now we take a closer look at Proposition 3.1. Part (i) provides a recursive
method to compute $P_{s}^{*}(\alpha)$ for all $s\in[n]$. For each $s\in[n]$,
by the expression of $\omega_{s}$, $P_{s}(i,\cdot)$ is a piecewise lower
semicontinuous linear-quadratic function whose domain is a closed interval,
relative to which $P_{s}(i,\cdot)$ has an expression of the form
$H(i)+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s}\|^{2}+(s-i)|\alpha|_{0}+\lambda_{1}$, while
$P_{s}^{*}(\cdot)=\min\\{P_{s}(0,\cdot),P_{s}(1,\cdot),\ldots,P_{s}(s-1,\cdot)\\}$.
Note that for each $i\in[0\\!:\\!s\\!-\\!1]$, the optimal solution to
$\min_{\alpha\in\mathbb{R}}P_{s}(i,\alpha)$ is easily obtained (in fact, all
the possible candidates of the global solutions are
$0,\frac{\sum_{j=i+1}^{s}z_{j}}{s-i},\max_{j\in[i+1:s]}\\{l_{j}\\},\min_{j\in[i+1:s]}\\{u_{j}\\}$),
so is
$\mathop{\arg\min}_{\alpha^{\prime}\in\mathbb{R}}P_{s}^{*}(\alpha^{\prime})$.
Part (ii) suggests a way to search for $i_{s}^{*}$ such that
$P_{s}^{*}(\alpha_{s}^{*})=P_{s}(i_{s}^{*},\alpha_{s}^{*})$ for each
$s\in[n]$. Obviously,
$P_{s}(i_{s}^{*},\alpha_{s}^{*})=\min_{i\in[0:s-1],\alpha\in\mathbb{R}}P_{s}(i,\alpha).$
This inspires us to propose Algorithm 1 for solving ${\rm
prox}_{\lambda_{1}\|\widehat{B}\cdot\|_{0}+\omega(\cdot)}(z)$, whose iterate
steps are described as follows.

Algorithm 1 (Computing ${\rm
prox}_{\lambda_{1}\|\widehat{B}\cdot\|_{0}+\omega(\cdot)}(z)$)

1\. Initialize: Compute
$P_{1}^{*}(\alpha)=\frac{1}{2}(z_{1}-\alpha)^{2}+\omega_{1}(\alpha)$ and set
$\mathcal{R}_{1}^{0}=\mathbb{R}.$

2\. For $s=2,\ldots,n\ $ do

3\.
$P_{s}^{*}(\alpha):=\min\\{P_{s-1}^{*}(\alpha),\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})+\lambda_{1}\\}+\frac{1}{2}(\alpha-
z_{s})^{2}+\omega_{s}(\alpha).$

4\. Compute $\mathcal{R}_{s}^{s-1}$ by (31).

5\. For $i=0,\ldots,s-2$ do

6\.
$\mathcal{R}_{s}^{i}=\mathcal{R}_{s-1}^{i}\cap(\mathcal{R}_{s}^{s-1})^{c}$.

7\. End

8\. End

9\. Set the current changepoint $s=n$.

10\. While $s>0$ do

11\. Find
$\alpha_{s}^{*}\in\mathop{\arg\min}_{\alpha\in\mathbb{R}}P_{s}^{*}(\alpha)$,
and $i_{s}^{*}=\left\\{i\ \lvert\
\alpha_{s}^{*}\in\mathcal{R}_{s}^{i}\right\\}$.

12\. $x^{*}_{i_{s}^{*}+1:s}=\alpha_{s}^{*}{\bf 1}$ and $s\leftarrow
i_{s}^{*}.$

13\. End

The main computation cost of Algorithm 1 comes from lines 3 and 6, in which
the number of pieces of the linear-quadratic functions involved in $P_{s}^{*}$
plays a crucial role. The following lemma gives a worst-case estimation for it
in the $s$-th iterate.

###### Lemma 3.4

Fix any $s\in[2:n]$. The function $P_{s}^{*}$ in line 3 of Algorithm 1 has at
most $O(s^{1+o(1)})$ linear-quadratic pieces.

Proof: Let $h_{i}(\alpha):=H(i)+\frac{1}{2}\|\alpha{\bf
1}-z_{i+1:s}\|^{2}+\lambda_{1}+(s-i)\lambda_{2}|\alpha|_{0}+\sum_{j=i+1}^{s}\delta_{[l_{j},u_{j}]}(\alpha)$
for $\alpha\in\mathbb{R}$ with $i\in[0\\!:\\!s\\!-\\!2]$. For each
$i\in[0\\!:\\!s\\!-\\!2]$, $h_{i}$ is a piecewise lower semicontinuous linear-
quadratic function whose domain is a closed interval, and every piece is
continuous on the closed interval except $\alpha=0$. Therefore, for each
$i\in[0\\!:\\!s\\!-\\!2]$,
$h_{i}=\min\big{\\{}h_{i,1},h_{i,2},h_{i,3}\big{\\}}$ with
$h_{i,1}(\alpha):=h_{i}(\alpha)-(s-i)\lambda_{2}|\alpha|_{0}+(s-i)\lambda_{2}+\delta_{(-\infty,0]}(\alpha),\,h_{i,2}(\alpha):=h_{i}(\alpha)+\delta_{\\{0\\}}(\alpha)$
and
$h_{i,3}(\alpha):=h_{i}(\alpha)-(s-i)\lambda_{2}|\alpha|_{0}+(s-i)\lambda_{2}+\delta_{[0,\infty)}(\alpha)$.
Obviously, $h_{i1},h_{i,2}$ and $h_{i,3}$ are piecewise linear-quadratic
functions with domain being a closed interval. In addition, write
$h_{s-1}(\alpha):=\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})+\frac{1}{2}(\alpha-
z_{s})^{2}+\lambda|\alpha|_{0}+\lambda_{1}+\delta_{[l_{s},u_{s}]}(\alpha)$ for
$\alpha\in\mathbb{R}$. Obviously, $h_{s-1}$ is a piecewise lower
semicontinuous linear-quadratic function whose domain is a closed interval.
From the above arguement, $h_{s-1}=\min\\{h_{s-1,1},h_{s-1,2},h_{s-1,3}\\}$
where each $h_{s-1,j}$ for $j=1,2,3$ is a piecewise linear function whose
domain is a closed interval. Combining the above discussion with line 3 of
Algorithm 1 and the definition of $P_{s-1}^{*}$, for any
$\alpha\in\mathbb{R}$,

$\displaystyle P_{s}^{*}(\alpha)\\!$
$\displaystyle=\min_{i\in[0:s-2]}\Big{\\{}P_{s-1}(i,\alpha)\\!+\\!\frac{1}{2}(\alpha-
z_{s})^{2}\\!+\\!\omega_{s}(\alpha),\min_{\alpha^{\prime}\in\mathbb{R}}P_{s-1}^{*}(\alpha^{\prime})\\!+\\!\frac{1}{2}(\alpha-
z_{s})^{2}\\!+\\!\omega_{s}(\alpha)\\!+\\!\lambda_{1}\Big{\\}}$
$\displaystyle=\Big{\\{}h_{0}(\alpha),h_{1}(\alpha),\ldots,h_{s-2}(\alpha),h_{s-1}(\alpha)\Big{\\}}=\min_{i\in[0:s-1],j\in[3]}\big{\\{}h_{i,j}(\alpha)\big{\\}}.$

Notice that any $h_{i,j}$ and $h_{i^{\prime},j^{\prime}}$ with $i\neq
i^{\prime}\in[0:s\\!-\\!1]$ or $j\neq j^{\prime}\in[3]$ crosses at most $2$
times. From [36, Theorem 2.5] the maximal number of linear-quadratic pieces
involved in $P_{s}^{*}$ is bounded by the maximal length of a $(3s,4)$
Davenport-Schinzel sequence, which by [14, Theorem 3] is
$3c_{1}s\exp(c_{2}\sqrt{\log 3s})$. Here, $c_{1},c_{2}$ are positive constants
independent of $s$. Thus, we conclude that the maximal number of linear-
quadratic pieces involved in $P_{s}^{*}$ is $O(s^{1+o(1)})$. The proof is
finished. $\Box$

By invoking Lemma 3.4, we are able to provide a worst-case estimation for the
complexity of Algorithm 1. Indeed, the main cost of Algorithm 1 consists in
lines 3 and 5-7. Since line 3 involves the computation cost proportional to
the pieces of $P_{s-1}^{*}$, from Lemma 3.4, it requires $O(s^{1+o(1)})$
operation. For each $i\in[0\\!:\\!s\\!-\\!1]$, from part (b) of Proposition
3.1 (ii), we know that $\mathcal{R}_{s}^{i}$ consists of at most
$O(s^{1+o(1)})$ intervals, which means that line 6 requires at most
$O(s^{1+o(1)})$ operations and then the complexity of lines 5-7 is
$O(s^{2+o(1)})$. Thus, the worst-case complexity of Algorithm 1 is
$\sum_{s=2}^{n}O(s^{2+o(1)})=O(n^{3+o(1)})$.

##  4 A hybrid of PG and inexact projected regularized Newton methods

In the hybrid frameworks owing to [40] and [5], the PG and Newton steps are
alternating. We now state the details of our algorithm, a hybrid of PG and
inexact projected regularized Newton methods (PGiPN), for solving problem (1),
where the introduction of the switch condition (8) is due to the consideration
that the PG step is more cost-effective than the Newton step when the iterates
are far from a stationary point. Let $x^{k}\in\Omega$ be the current iterate.
It is noted that the PG step is always executed and if condition (8) is met,
we need to solve (11), which involves constructing $G_{k}$ to satisfy (12).
Such $G^{k}$ can be easily achieved in the following situations.

For some generalized linear models, $f$ can be expressed as $f(x)=h(Ax-b)$ for
some $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$ and twice continuously
differentiable, separable $h$. For this case, $\nabla^{2}h$ is a diagonal
matrix, and $\nabla^{2}\\!f(x)=A^{\top}\nabla^{2}h(Ax-b)A$. Since
$\nabla^{2}\\!f(x^{k})$ is not necessarily positive definite, following the
method in [26], we construct $G_{k}:=G_{k}^{1}$, where

$G_{k}^{1}:=\nabla^{2}f(x^{k})\\!+b_{2}[-\lambda_{\min}(\nabla^{2}\\!h(Ax^{k}-b))]_{+}A^{\top}A+b_{1}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}I$
(34)

with $b_{2}\geq 1$. However, for highly nonconvex $h$,
$[-\lambda_{\min}(\nabla^{2}\\!h(Ax^{k}-b))]_{+}$ is large, for which
$G_{k}^{1}$ is a poor approximation to $\nabla^{2}f(x^{k})$. To avoid this
drawback and simultaneously make $G_{k}$ positive definite, we consider
$G_{k}:=G_{k}^{2}$, where

$G_{k}^{2}:=A^{\top}[\nabla^{2}\\!h(Ax^{k}-b)]_{+}A+b_{1}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}I.$
(35)

For the case where $\nabla^{2}\\!f(x^{k})\succeq 0$, $G_{k}^{1}=G_{k}^{2}$. If
$\nabla^{2}\\!f(x^{k})\not\succeq 0$, it is immediate to see that
$\|G_{k}^{1}-\nabla^{2}\\!f(x^{k})\|_{2}\geq\|G_{k}^{2}-\nabla^{2}\\!f(x^{k})\|_{2}$,
which means that $G_{k}^{2}$ is a better approximation to $\nabla^{2}f(x^{k})$
than $G_{k}^{1}$. On the other hand, for those $f$’s not owning a separable
structure, we form $G_{k}:=G_{k}^{3}$ as in [42] and [45], where

$G_{k}^{3}:=\nabla^{2}f(x^{k})\\!+\big{(}b_{2}[-\lambda_{\min}(\nabla^{2}f(x^{k}))]_{+}+b_{1}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}\big{)}I.$
(36)

It is not hard to check that for $i=1,2,3$, $G_{k}^{i}$ meets the requirement
in (12). We remark here that the sequel convergence analysis holds for all
three $G_{k}^{i}$, and we write them by $G_{k}$ for simplicity.

The iterates of PGiPN are described as follows.

Algorithm 2 (a hybrid of PG and inexact projected regularized Newton methods)

Initialization: Choose $\epsilon\geq 0$ and parameters
$\mu_{\max}>\mu_{\min}>0,\,\tau>1,\,\alpha>0,b_{1}>0,b_{2}\geq 1$,
$\varrho\in(0,\frac{1}{2}),\,\sigma\in(0,\frac{1}{2}),\,\varsigma\in(\sigma,1]$
and $\beta\in(0,1)$. Choose an initial $x^{0}\in\Omega$ and let $k:=0$.

PG Step:

(1a) Select $\mu_{k}\in[\mu_{\min},\mu_{\max}]$. Let $m_{k}$ be the smallest
nonnegative integer $m$ such that

$F(\overline{x}^{k})\leq
F(x^{k})-\frac{\alpha}{2}\|x^{k}\\!-\\!\overline{x}^{k}\|^{2}\ \ {\rm with}\ \
\overline{x}^{k}\in{\rm
prox}_{(\mu_{k}\tau^{m})^{-1}g}(x^{k}\\!-\\!(\mu_{k}\tau^{m})^{-1}\nabla
f(x^{k})).$ (37)

(1b) Let $\overline{\mu}_{k}=\mu_{k}\tau^{m_{k}}$. If
$\overline{\mu}_{k}\|x^{k}-\overline{x}^{k}\|\leq\epsilon,$ stop and output
$x^{k}$; otherwise, go to step (1c).

(1c) If condition (8) holds, go to Newton step; otherwise, let
$x^{k+1}=\overline{x}^{k}$. Set $k\leftarrow k+1$ and return to step (1a).

Newton step:

(2a) Seek an inexact solution $y^{k}$ of (11) satisfying (13)-(14).

(2b) Set $d^{k}:=y^{k}-x^{k}$. Let $t_{k}$ be the smallest nonnegative integer
$t$ such that

$f(x^{k}+\beta^{t}d^{k})\leq f(x^{k})+\varrho\beta^{t}\langle\nabla
f(x^{k}),d^{k}\rangle.$ (38)

(2c) Let $\alpha_{k}=\beta^{t_{k}}$ with
$x^{k+1}=x^{k}\\!+\\!\alpha_{k}d^{k}$. Set $k\leftarrow k+1$ and return to PG
step.

###### Remark 4.1

(i) Our PGiPN benefits from the PG step in two aspects. First, the
incorporation of the PG step can guarantee that the sequence generated by
PGiPN remains in a right position for convergence. Second, the PG step helps
to identify adaptively the subspace used in the Newton step, and as will be
shown in Proposition 5.1, switch condition (8) always holds and the supports
of $\\{Bx^{k}\\}_{k\in\mathbb{N}}$ and $\\{x^{k}\\}_{k\in\mathbb{N}}$ keep
unchanged when $k$ is sufficiently large, so that Algorithm 2 will reduce to
an inexact projected regularized Newton method for solving (10) with
$\Pi_{k}\equiv\Pi_{*}$. In this sense, the PG step plays a crucial role in
transforming the original challenging problem (1) into a problem that can be
efficiently solved by the inexact projected regularized Newton method.

(ii) When $x^{k}$ enters the Newton step, from the inexact criterion (13) and
the expression of $\Theta_{k}$,
$0\geq\Theta_{k}(x^{k}\\!+\\!d^{k})-\Theta_{k}(x^{k})=\langle\nabla
f(x^{k}),d^{k}\rangle+\frac{1}{2}\langle d^{k},G_{k}d^{k}\rangle$, and then

$\langle\nabla f(x^{k}),d^{k}\rangle\leq-\frac{1}{2}\langle
d^{k},G_{k}d^{k}\rangle\leq-\frac{b_{1}}{2}\|\overline{\mu}_{k}(x^{k}-\overline{x}^{k})\|^{\sigma}\|d^{k}\|^{2}<0,$
(39)

where the second inequality is due to (12). In addition, the inexact criterion
(13) implies that $y^{k}\in\Pi_{k}$, which along with $x^{k}\in\Pi_{k}$ and
the convexity of $\Pi_{k}$ yields that $x^{k}+\alpha d^{k}\in\Pi_{k}$ for any
$\alpha\in(0,1]$. By the definition of $\Pi_{k}$, ${\rm supp}(B(x^{k}+\alpha
d^{k}))\subset{\rm supp}(Bx^{k})$ and ${\rm supp}(x^{k}+\alpha
d^{k})\subset{\rm supp}(x^{k})$, so $g(x^{k}+\alpha d^{k})\leq g(x^{k})$ for
any $\alpha\in(0,1]$. This together with (39) shows that the iterate along the
direction $d^{k}$ will reduce the value of $F$ at $x^{k}$.

(iii) When $\epsilon=0$, by Definition 2.1 the output $x^{k}$ of Algorithm 2
is an $L$-stationary point of (1), which is also a stationary point of problem
(10) from Proposition 2.1 and Lemma 2.1 (i). Let
$r_{k}\\!:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be the KKT residual mapping of (10)
defined by

$r_{k}(x):=\overline{\mu}_{k}[x-{\rm
proj}_{\Pi_{k}}(x-\overline{\mu}_{k}^{-1}\nabla f(x))].$ (40)

It is not difficult to verify that when $x^{k}$ satisfies condition (8), the
following relation holds

$r_{k}(x^{k})=\overline{\mu}_{k}(x^{k}-\overline{x}^{k}),$ (41)

for which it suffices to argue that $\overline{x}^{k}\\!={\rm
proj}_{\Pi_{k}}(x^{k}-\overline{\mu}_{k}^{-1}\nabla f(x^{k}))$. Indeed, if
not, there exists $\overline{z}^{k}\in\Pi_{k}$ such that
$\widetilde{h}_{k}(\overline{z}^{k})<\widetilde{h}_{k}(\overline{x}^{k})$,
where
$\widetilde{h}_{k}(x):=\frac{\overline{\mu}_{k}}{2}\|x-(x^{k}-\overline{\mu}_{k}^{-1}\nabla
f(x^{k}))\|^{2}$. Since $\overline{z}^{k}\in\Pi_{k}$, we have ${\rm
supp}(B\overline{z}^{k})\subset{\rm supp}(B\overline{x}^{k})$ and ${\rm
supp}(\overline{z}^{k})\subset{\rm supp}(\overline{x}^{k})$, which implies
that $g(\overline{z}^{k})\leq g(\overline{x}^{k})$ and then
$\widetilde{h}_{k}(\overline{z}^{k})+g(\overline{z}^{k})<\widetilde{h}_{k}(\overline{x}^{k})+g(\overline{x}^{k})$,
a contradiction to $\overline{x}^{k}\in{\rm
prox}_{\overline{\mu}_{k}^{-1}g}(x^{k}-\overline{\mu}_{k}^{-1}\nabla
f(x^{k}))$.

(iv) By using (16) and the descent lemma [8, Proposition A.24], the line
search in step (1a) must stop after a finite number of backtrackings. In fact,
the line search in step (1a) is satisfied when $\mu_{k}\tau^{m}\geq
L_{1}+\alpha$, which implies that
$\overline{\mu}_{k}<\widetilde{\mu}:=\tau(L_{1}\\!+\alpha)$ for each
$k\in\mathbb{N}$.

By Remark 4.1 (iv), to show that Algorithm 2 is well defined, we only need to
argue that the Newton steps in Algorithm 2 are well defined, which is implied
by the following lemma.

###### Lemma 4.1

For each $k\in\mathbb{N}$, define the KKT residual mapping
$R_{k}\\!:\mathbb{R}^{n}\to\mathbb{R}^{n}$ of (11) by

$R_{k}(y):=\overline{\mu}_{k}[y-{\rm
proj}_{\Pi_{k}}(y-\overline{\mu}_{k}^{-1}(G_{k}(y-x^{k})+\nabla\\!f(x^{k})))].$

Then, for those $x^{k}$’s satisfying (8), the following statements are true.

(i)

For any $y$ close enough to the optimal solution of (11),
$y-\overline{\mu}_{k}^{-1}R_{k}(y)$ satisfies (13)-(14).

(ii)

The line search step in (38) is well defined, and
$\alpha_{k}\geq\min\big{\\{}1,\frac{(1-\varrho)b_{1}\beta}{L_{1}}\|\overline{\mu}_{k}(x^{k}-\overline{x}^{k})\|^{\sigma}\big{\\}}$.

(iii)

The inexact criterion (14) implies that
$\|R_{k}(y^{k})\|\leq\frac{1}{2}\min\left\\{\|r_{k}(x^{k})\|,\|r_{k}(x^{k})\|^{1+\varsigma}\right\\}$.

Proof: Pick any $x^{k}$ satisfying (8). We proceed the proof of parts
(i)-(iii) as follows.

(i) Let $\widehat{y}^{k}$ be the unique optimal solution to (11). Then
$\widehat{y}^{k}\neq x^{k}$ (if not, $x^{k}$ is the optimal solution of (11)
and $0=R_{k}(x^{k})=r_{k}(x^{k})$, which by (41) means that
$x^{k}=\overline{x}^{k}$ and Algorithm 2 stops at $x^{k}$). By the optimality
condition of (11),
$-\nabla\\!f(x^{k})-G_{k}(\widehat{y}^{k}-x^{k})\in\mathcal{N}_{\Pi_{k}}(\widehat{y}^{k})$,
which by the convexity of $\Pi_{k}$ and $x^{k}\in\Pi_{k}$ implies that
$\langle\nabla\\!f(x^{k})+G_{k}(\widehat{y}^{k}-x^{k}),\widehat{y}^{k}-x^{k}\rangle\leq
0$. Along with the expression of $\Theta_{k}$, we have
$\Theta_{k}(\widehat{y}^{k})-\Theta_{k}(x^{k})\leq-\frac{1}{2}\langle\widehat{y}^{k}-x^{k},G_{k}(\widehat{y}^{k}-x^{k})\rangle<0$.
Since $\Theta_{k}$ is continuous relative to $\Pi_{k}$, for any $z\in\Pi_{k}$
sufficiently close to $\widehat{y}^{k}$, $\Theta_{k}(z)-\Theta_{k}(x^{k})\leq
0$. From $R_{k}(\widehat{y})=0$ and the continuity of $R_{k}$,
$y-\overline{\mu}_{k}^{-1}R_{k}(y)$ is close to $\widehat{y}$ when $y$
sufficiently close to $\widehat{y}$, which together with
$y-\overline{\mu}_{k}^{-1}R_{k}(y)\in\Pi_{k}$ implies that
$y-\overline{\mu}_{k}^{-1}R_{k}(y)$ satisfies the criterion (13) when $y$ is
sufficiently close to $\widehat{y}$. In addition, from the expression of
$R_{k}$, for any $y\in\mathbb{R}^{n}$,

$0\in
G_{k}(y-x^{k})+\nabla\\!f(x^{k})-R_{k}(y)+\mathcal{N}_{\Pi_{k}}(y-\overline{\mu}_{k}^{-1}R_{k}(y)),$

which by the expression of $\Theta_{k}$ implies that
$\overline{\mu}_{k}^{-1}G_{k}R_{k}(y)+R_{k}(y)\in\partial\Theta_{k}(y-\overline{\mu}_{k}^{-1}R_{k}(y))$.
Hence, ${\rm
dist}(0,\partial\Theta_{k}(y-\overline{\mu}_{k}^{-1}R_{k}(y)))\leq\|\overline{\mu}_{k}^{-1}G_{k}R_{k}(y)+R_{k}(y)\|$.
Noting that $R_{k}(\widehat{y}^{k})=0$, we have
$\|\overline{\mu}_{k}^{-1}G_{k}R_{k}(\widehat{y}^{k})+R_{k}(\widehat{y}^{k})\|=0<\frac{\min\\{\overline{\mu}_{k}^{-1},1\\}}{2}\min\left\\{\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|,\|\overline{\mu}_{k}(x^{k}\\!-\overline{x}^{k})\|^{1+\varsigma}\right\\}$.
From the continuity of the function
$y\mapsto\|\overline{\mu}_{k}^{-1}G_{k}R_{k}(y)+R_{k}(y)\|$, we conclude that
for any $y$ sufficiently close to $\widehat{y}^{k}$,
$y-\overline{\mu}_{k}^{-1}R_{k}(y)$ satisfies the inexact criterion (14).

(ii) By (16) and the descent lemma [8, Proposition A.24], for any
$\alpha\in(0,1]$,

$\displaystyle f(x^{k}\\!+\\!\alpha
d^{k})-f(x^{k})-\varrho\alpha\langle\nabla\\!f(x^{k}),d^{k}\rangle$
$\displaystyle\leq(1\\!-\\!\varrho)\alpha\langle\nabla\\!f(x^{k}),d^{k}\rangle+\frac{L_{1}\alpha^{2}}{2}\|d^{k}\|^{2}$
$\displaystyle\leq-\frac{(1\\!-\\!\varrho)\alpha
b_{1}}{2}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}\|d^{k}\|^{2}\\!+\\!\frac{L_{1}\alpha^{2}}{2}\|d^{k}\|^{2}$
$\displaystyle=\Big{(}-\frac{(1\\!-\\!\varrho)b_{1}}{2}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}\\!+\\!\frac{L_{1}\alpha}{2}\Big{)}\alpha\|d^{k}\|^{2},$

where the second inequality uses (39). Therefore, when
$\alpha\leq\min\big{\\{}1,\frac{(1-\varrho)b_{1}}{L_{1}}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}\big{\\}}$,
the line search in (38) holds, which implies that
$\alpha_{k}\geq\min\big{\\{}1,\frac{(1-\varrho)b_{1}\beta}{L_{1}}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}\big{\\}}$.

(iii) Let $\zeta^{k}\in\partial\Theta_{k}(y^{k})$ be such that
$\|\zeta^{k}\|={\rm dist}(0,\partial\Theta_{k}(y^{k}))$. From
$\zeta^{k}\in\partial\Theta_{k}(y^{k})$ and the expression of $\Theta_{k}$, we
have $y^{k}={\rm proj}_{\Pi_{k}}(y^{k}+\zeta^{k}-(G_{k}(y^{k}-x^{k})+\nabla
f(x^{k})))$. Along with the nonexpansiveness of ${\rm proj}_{\Pi_{k}}$,
$\|y^{k}\\!-\\!{\rm
proj}_{\Pi_{k}}(y^{k}\\!-\\!(G_{k}(y^{k}\\!-\\!x^{k})+\nabla
f(x^{k})))\|\leq\|\zeta^{k}\|$. Consequently,

$\displaystyle{\rm dist}(0,\partial\Theta_{k}(y^{k}))$
$\displaystyle\geq\|y^{k}\\!-\\!{\rm
proj}_{\Pi_{k}}(y^{k}\\!-\\!(G_{k}(y^{k}\\!-\\!x^{k})+\nabla
f(x^{k})))\|\geq\min\\{\overline{\mu}_{k}^{-1},1\\}\|R_{k}(y^{k})\|,$

where the second inequality follows by Lemma 4 of [37] and the expression of
$R_{k}$. Combining the last inequality with (14) and (41) leads to the desired
inequality. $\Box$

When $\overline{\mu}_{k}=1$, the condition that
$\|R_{k}(y^{k})\|\leq\frac{1}{2}\min\left\\{\|r_{k}(x^{k})\|,\|r_{k}(x^{k})\|^{1+\varsigma}\right\\}$
is a special case of the first inexact condition in [46, Equa (6a)] or the
inexact condition in [30, Equa (14)], which by Lemma 4.1 (iii) shows that
criterion (14) with $\overline{\mu}_{k}=1$ is stronger than those ones.

To analyze the convergence of Algorithm 2 with $\epsilon=0$, henceforth we
assume $x^{k}\neq\overline{x}^{k}$ for all $k$ (if not, Algorithm 2 will
produce an $L$-stationary point within finite number of steps, and its
convergence holds automatically). From the iterate steps of Algorithm 2, we
see that the sequence $\\{x^{k}\\}_{k\in\mathbb{N}}$ consists of two parts,
$\\{x^{k}\\}_{k\in\mathcal{K}_{1}}$ and $\\{x^{k}\\}_{k\in\mathcal{K}_{2}}$,
where

$\mathcal{K}_{1}\\!:=\\!\mathbb{N}\backslash\mathcal{K}_{2}\ \ {\rm with}\ \
\mathcal{K}_{2}\\!:=\\!\big{\\{}k\in\mathbb{N}\ \lvert\ {\rm
supp}(Bx^{k})\\!={\rm supp}(B\overline{x}^{k}),\ {\rm supp}(x^{k})\\!={\rm
supp}(\overline{x}^{k})\big{\\}}.$

Obviously, $\mathcal{K}_{1}$ consists of those $k$’s with $x^{k+1}$ from the
PG step, while $\mathcal{K}_{2}$ consists of those $k$’s with $x^{k+1}$ from
the Newton step.

To close this section, we provide some properties of the sequences
$\\{x^{k}\\}_{k\in\mathbb{N}}$ and $\\{\overline{x}^{k}\\}_{k\in\mathbb{N}}$.

###### Proposition 4.1

The following assertions are true.

(i)

The sequence $\\{F(x^{k})\\}_{k\in\mathbb{N}}$ is descent and convergent.

(ii)

There exists $\nu>0$ such that $\lvert B\overline{x}^{k}\rvert_{\min}\geq\nu$
and $\lvert\overline{x}^{k}\rvert_{\min}\geq\nu$ for all $k\in\mathbb{N}$.

(iii)

There exist $c_{1},c_{2}>0$ such that $c_{1}\|r_{k}(x^{k})\|\leq\|d^{k}\|\leq
c_{2}\|r_{k}(x^{k})\|^{1-\sigma}$ for all $k\in\mathcal{K}_{2}$.

Proof: (i) For each $k\in\mathbb{N}$, when $k\in\mathcal{K}_{1}$, by the line
search in step (1a), $F(x^{k+1})<F(x^{k})$, and when $k\in\mathcal{K}_{2}$,
from (38) and (39), it follows that $f(x^{k+1})<f(x^{k})$, which along with
$g(x^{k+1})\leq g(x^{k})$ by Remark 4.1 (ii) implies that
$F(x^{k+1})<F(x^{k})$. Hence, $\\{F(x^{k})\\}_{k\in\mathbb{N}}$ is a descent
sequence. Recall that $F$ is lower bounded on $\Omega$, so
$\\{F(x^{k})\\}_{k\in\mathbb{N}}$ is convergent.

(ii) By the definition of $\overline{\mu}_{k}$ and Remark 4.1 (iv),
$\overline{\mu}_{k}\in[\mu_{\min},\widetilde{\mu})$ for all $k\in\mathbb{N}$.
Note that $\\{x^{k}\\}_{k\in\mathbb{N}}\subset\Omega$, so the sequence
$\\{x^{k}\\!-\overline{\mu}_{k}^{-1}\nabla f(x^{k})\\}_{k\in\mathbb{N}}$ is
bounded and is contained in a compact set, says, $\Xi$. By invoking Lemma 3.2
with such $\Xi$ and
$\underline{\mu}=\mu_{\min},\overline{\mu}=\widetilde{\mu}$, there exists
$\nu>0$ (depending on $\Xi,\ \mu_{\min}$ and $\widetilde{\mu}$) such that
$\lvert[B;I]\overline{x}^{k}\rvert_{\min}>\nu$. The desired result then
follows by noting that $\lvert
B\overline{x}^{k}\rvert_{\min}\geq\lvert[B;I]\overline{x}^{k}\rvert_{\min}$
and
$\lvert\overline{x}^{k}\rvert_{\min}\geq\lvert[B;I]\overline{x}^{k}\rvert_{\min}$.

(iii) From the definition of $G_{k}$, the continuity of $\nabla^{2}\\!f$,
$\\{x^{k},\overline{x}^{k}\\}_{k\in\mathbb{N}}\subset\Omega$ and Remark 4.1
(iv), there exists $\overline{c}>0$ such that

$\|G_{k}\|_{2}\leq\overline{c}\ {\rm for\ all}\ k\in\mathcal{K}_{2}.$ (42)

Fix any $k\in\mathcal{K}_{2}$. By Lemma 4.1 (iii),
$\|R_{k}(y^{k})\|\leq\frac{1}{2}\|r_{k}(x^{k})\|$. Then, it holds that

$\displaystyle\frac{1}{2}\|r_{k}(x^{k})\|\leq\|r_{k}(x^{k})\|-\|R_{k}(y^{k})\|\leq\|r_{k}(x^{k})-R_{k}(y^{k})\|$
$\displaystyle=\overline{\mu}_{k}\|x^{k}-{\rm
proj}_{\Pi_{k}}(x^{k}-\overline{\mu}_{k}^{-1}\nabla f(x^{k}))-y^{k}+{\rm
proj}_{\Pi_{k}}(y^{k}-\overline{\mu}_{k}^{-1}(G_{k}(y^{k}-x^{k})+\nabla\\!f(x^{k})))\|$
$\displaystyle\leq(2\overline{\mu}_{k}+\|G_{k}\|_{2})\|y^{k}-x^{k}\|\leq(2\widetilde{\mu}+\overline{c})\|d^{k}\|,$

where the third inequality is using the nonexpansiveness of ${\rm
proj}_{\Pi_{k}}$, and the last one is due to (42) and $d^{k}=y^{k}-x^{k}$.
Therefore, $c_{1}\|r_{k}(x^{k})\|\leq\|d^{k}\|$ with
$c_{1}:=1/(4\widetilde{\mu}\\!+2\overline{c})$. For the second inequality, it
follows from the definitions of $r_{k}(\cdot)$ and $R_{k}(\cdot)$ that
$R_{k}(y^{k})-\nabla
f(x^{k})-G_{k}d^{k}\in\mathcal{N}_{\Pi_{k}}(y^{k}-\overline{\mu}_{k}^{-1}R_{k}(y^{k}))$
and $r_{k}(x^{k})-\nabla
f(x^{k})\in\mathcal{N}_{\Pi_{k}}(x^{k}-\overline{\mu}_{k}^{-1}r_{k}(x^{k}))$,
which together with the monotonicity of the set-valued mapping
$\mathcal{N}_{\Pi_{k}}(\cdot)$ implies that

$\displaystyle\langle d^{k},G_{k}d^{k}\rangle$ $\displaystyle\leq\langle
R_{k}(y^{k})\\!-\\!r_{k}(x^{k}),d^{k}\rangle-\overline{\mu}_{k}^{-1}\|R_{k}(y^{k})\\!-\\!r_{k}(x^{k})\|^{2}-\overline{\mu}_{k}^{-1}\langle
G_{k}d^{k},-R_{k}(y^{k})+r_{k}(x^{k})\rangle$
$\displaystyle\leq\langle(I+\overline{\mu}_{k}^{-1}G_{k})d^{k},R_{k}(y^{k})-r_{k}(x^{k})\rangle.$

Combining this inequality with equations (12), (41) and Lemma 4.1 (iii) leads
to

$\displaystyle b_{1}\|r_{k}(x^{k})\|^{\sigma}\|d^{k}\|^{2}$
$\displaystyle\leq(1+\overline{\mu}_{k}^{-1}\|G_{k}\|_{2})(\|R_{k}(y^{k})\|+\|r_{k}(x^{k})\|)\|d^{k}\|$
(43)
$\displaystyle\leq(3/2)(1+\overline{\mu}_{k}^{-1}\|G_{k}\|_{2})\|r_{k}(x^{k})\|\|d^{k}\|,$

which along with (42) and $\mu_{k}\geq\mu_{\min}$ implies that
$\|d^{k}\|\leq\frac{3}{2}(1+\mu_{\min}^{-1}\overline{c})b_{1}^{-1}\|r_{k}(x^{k})\|^{1-\sigma}.$
Then, $\|d^{k}\|\leq c_{2}\|r_{k}(x^{k})\|^{1-\sigma}$ holds with
$c_{2}:=\frac{3}{2}(1+\mu_{\min}^{-1}\overline{c})b_{1}^{-1}$. The proof is
completed. $\Box$

##  5 Convergence Analysis

Before analyzing the convergence of Algorithm 2, we show that Algorithm 2
finally reduces to an inexact projected regularized Newton method for seeking
a stationary point of a problem to minimize a smooth function over a
polyhedral set. This requires the following lemma.

###### Lemma 5.1

For the sequences $\\{x^{k}\\}_{k\in\mathbb{N}}$ and
$\\{\overline{x}^{k}\\}_{k\in\mathbb{N}}$ generated by Algorithm 2, the
following assertions are true.

(i)

There exists a constant $\gamma>0$ such that for each $k\in\mathbb{N}$,

$F(x^{k+1})-F(x^{k})\leq\left\\{\begin{array}[]{ll}-\gamma\|x^{k}-\overline{x}^{k}\|^{2}&{\rm
if}\ k\in\mathcal{K}_{1},\\\ -\gamma\|x^{k}-\overline{x}^{k}\|^{2+\sigma}&{\rm
if}\ k\in\mathcal{K}_{2},\ \alpha_{k}=1,\\\
-\gamma\|x^{k}-\overline{x}^{k}\|^{2+2\sigma}&{\rm if}\ k\in\mathcal{K}_{2},\
\alpha_{k}\neq 1.\end{array}\right.$ (44)

(ii)

$\lim_{k\rightarrow\infty}\|x^{k}-\overline{x}^{k}\|=0$ and
$\lim_{\mathcal{K}_{2}\ni k\rightarrow\infty}\|d^{k}\|=0$.

(iii)

The accumulation point set of $\\{x^{k}\\}_{k\in\mathbb{N}}$, denoted by
$\Gamma(x^{0})$, is nonempty and compact, and every element of $\Gamma(x^{0})$
is an $L$-stationary point of problem (1).

Proof: (i) Fix any $k\in\mathcal{K}_{2}$. From inequalities (38) and (39),

$\displaystyle f(x^{k+1})-f(x^{k})$ $\displaystyle\leq-\frac{\varrho
b_{1}\alpha_{k}}{2}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{\sigma}\|d^{k}\|^{2}\leq-\frac{\varrho
c_{1}^{2}b_{1}\alpha_{k}}{2}\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{2+\sigma}$
(45) $\displaystyle\leq-\frac{\varrho
c_{1}^{2}b_{1}\alpha_{k}\mu_{\min}^{2+\sigma}}{2}\|x^{k}\\!-\\!\overline{x}^{k}\|^{2+\sigma},$

where the second inequality is using Proposition 4.1 (iii) and equality (41).
By Remark 4.1 (ii), we have $g(x^{k+1})\leq g(x^{k})$, so that
$F(x^{k+1})-F(x^{k})\leq f(x^{k+1})-f(x^{k})$, which along with the last
equation yields that $F(x^{k+1})\\!-F(x^{k})\leq-\frac{\varrho
c_{1}^{2}b_{1}\alpha_{k}\mu_{\min}^{2+\sigma}}{2}\|x^{k}-\overline{x}^{k}\|^{2+\sigma}$.
Take $\gamma\\!:=\min\big{\\{}\frac{\alpha}{2},\frac{\varrho
c_{1}^{2}b_{1}\mu_{\min}^{2+\sigma}}{2},\frac{\beta(1-\varrho)\varrho
c_{1}^{2}b_{1}^{2}\mu_{\min}^{2+2\sigma}}{2L_{1}}\big{\\}}$. The desired
result then follows by using Lemma 4.1 (ii) and recalling that
$F(x^{k+1})-F(x^{k})\leq\frac{\alpha}{2}\|x^{k}-\overline{x}^{k}\|^{2}$ for
$k\in\mathcal{K}_{1}$.

(ii) Let $\widetilde{\mathcal{K}}_{2}:=\\{k\in\mathcal{K}_{2}\ |\
\alpha_{k}=1\\}$. Doing summation for inequality (44) from $i=1$ to any
$j\in\mathbb{N}$ yields that

$\displaystyle\sum_{i\in\mathcal{K}_{1}\cap[j]}\gamma\|x^{i}-\overline{x}^{i}\|^{2}$
$\displaystyle+\sum_{i\in\widetilde{\mathcal{K}}_{2}\cap[j]}\gamma\|x^{i}-\overline{x}^{i}\|^{2+\sigma}+\sum_{i\in(\mathcal{K}_{2}\backslash\widetilde{\mathcal{K}}_{2})\cap[j]}\gamma\|x^{i}-\overline{x}^{i}\|^{2+2\sigma}$
$\displaystyle\leq\sum_{i=1}^{j}\big{[}F(x^{i})-F(x^{i+1})\big{]}=F(x^{1})-F(x^{j+1}),$

which by the lower boundedness of $F$ on the set $\Omega$ implies that

$\sum_{i\in\mathcal{K}_{1}}\|x^{i}-\overline{x}^{i}\|^{2}+\sum_{i\in\widetilde{\mathcal{K}}_{2}}\gamma\|x^{i}-\overline{x}^{i}\|^{2+\sigma}+\sum_{i\in\mathcal{K}_{2}\backslash\widetilde{\mathcal{K}}_{2}}\gamma\|x^{i}-\overline{x}^{i}\|^{2+2\sigma}<\infty.$

Thus, we obtain $\lim_{k\rightarrow\infty}\|x^{k}-\overline{x}^{k}\|=0$.
Together with Proposition 4.1 (iii), (41) and Remark 4.1 (iv), it follows that
$\lim_{\mathcal{K}_{2}\ni k\rightarrow\infty}\|d^{k}\|=0$.

(iii) Recall that $\\{x^{k}\\}_{k\in\mathbb{N}}\subset\Omega$, so its
accumulation point set $\Gamma(x^{0})$ is nonempty. Pick any
$x^{*}\in\Gamma(x^{0})$. Then, there exists an index set $K\subset\mathbb{N}$
such that $\lim_{K\ni k\to\infty}x^{k}=x^{*}$. From part (ii), $\lim_{K\ni
k\to\infty}\overline{x}^{k}=x^{*}$. For each $k\in K$,
$\overline{x}^{k}\in{\rm
prox}_{\overline{\mu}_{k}^{-1}g}\big{(}x^{k}-\overline{\mu}_{k}^{-1}\nabla
f(x^{k})\big{)}$ with $\overline{\mu}_{k}\in[\mu_{\rm min},\widetilde{\mu})$
by step (1a) and Remark 4.1 (iv). We assume that $\lim_{K\ni
k\to\infty}\overline{\mu}_{k}\to\overline{\mu}_{*}\in[\mu_{\rm
min},\widetilde{\mu}]$ (if necessary by taking a subsequence). Define the
function

$h(z,x,\mu)\\!:=\\!\left\\{\begin{array}[]{cl}\frac{\mu}{2}\big{\|}z-(x\\!-\\!\mu^{-1}\nabla\\!f(x))\big{\|}^{2}+g(z)&{\rm
if}\ (z,x,\mu)\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times[\mu_{\rm
min},\widetilde{\mu}],\\\ \infty&{\rm otherwise},\end{array}\right.$

and write
$\mathcal{P}(x,\mu):=\mathop{\arg\min}_{z\in\mathbb{R}^{n}}h(z,x,\mu)$ for
$(x,\mu)\in\mathbb{R}^{n}\times\mathbb{R}$. Note that
$h\\!:\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}\to\overline{\mathbb{R}}$
is a proper and lower semicontinuous function and is level-bounded in $z$
locally uniformly in $(x,\mu)$. In addition, by [34, Theorem 1.25], the
function $\widehat{h}(x,\mu):=\inf_{z\in\mathbb{R}^{n}}h(z,x,\mu)$ is finite
and continuous on $\mathbb{R}\times[\mu_{\rm min},\widetilde{\mu}]$. From [34,
Example 5.22], the multifunction
$\mathcal{P}\\!:\mathbb{R}^{n}\times\mathbb{R}\rightrightarrows\mathbb{R}^{n}$
is outer semicontinuous relative to $\mathbb{R}^{n}\times[\mu_{\rm
min},\widetilde{\mu}]$. Note that
$\overline{x}^{k}\in\mathcal{P}(x^{k},\overline{\mu}_{k})$ for each $k\in K$.
Then, $x^{*}\in{\rm
prox}_{\overline{\mu}_{*}^{-1}g}\big{(}x^{*}\\!-\\!\overline{\mu}_{*}^{-1}\nabla
f(x^{*})\big{)}$, so $x^{*}$ is an $L$-stationary point of (1). $\Box$

Next we use Lemma 5.1 (ii) to show that, after a finite number of iterations,
the switch condition in (8) always holds and the Newton step is executed. To
this end, define

$T_{k}\\!:={\rm supp}(Bx^{k}),\ \overline{T}_{k}\\!:={\rm
supp}(B\overline{x}^{k}),\ S_{k}\\!:={\rm supp}(x^{k})\ \ {\rm and}\ \
\overline{S}_{k}\\!:={\rm supp}(\overline{x}^{k}).$ (46)

###### Proposition 5.1

For the index sets defined in (46), there exist index sets
$T\subset[p],S\subset[n]$ and an index $\overline{k}\in\mathbb{N}$ such that
for all $k>\overline{k}$, $T_{k}=\overline{T}_{k}=T$ and
$S_{k}=\overline{S}_{k}=S$, which means that $k\in\mathcal{K}_{2}$ for all
$k>\overline{k}$. Moreover, for each $x^{*}\in\Gamma(x^{0})$, ${\rm
supp}(Bx^{*})=T$ and ${\rm supp}(x^{*})=S$, where $\Gamma(x^{0})$ is defined
in Lemma 5.1 (iii).

Proof: We complete the proof of the conclusion via the following three claims:

Claim 1: There exists $\overline{k}\in\mathbb{N}$ such that for
$k>\overline{k}$, $\lvert Bx^{k}\rvert_{\min}\geq\frac{\nu}{2}$, where $\nu$
is the same as the one in Proposition 4.1 (ii). Indeed, for each
$k\\!-\\!1\in\mathcal{K}_{1}$, $x^{k}=\overline{x}^{k-1}$, and $\lvert
Bx^{k}\rvert_{\min}=\lvert
B\overline{x}^{k-1}\rvert_{\min}\geq\nu>\frac{\nu}{2}$ follows by Proposition
4.1 (ii). Hence, it suffices to consider that $k\\!-\\!1\in\mathcal{K}_{2}$.
By Lemma 5.1 (ii), there exists $\overline{k}\in\mathbb{N}$ such that for all
$k>\overline{k}$, $\|x^{k-1}-\overline{x}^{k-1}\|<\frac{\nu}{4\|B\|_{2}}$ and
$\|d^{k-1}\|<\frac{\nu}{4\|B\|_{2}}$, which implies that for
$\mathcal{K}_{2}\ni k-1>\overline{k}-1$,
$\|Bx^{k-1}\\!-\\!B\overline{x}^{k-1}\|<\frac{\nu}{4}$ and
$\|Bd^{k-1}\|<\frac{\nu}{4}$. For each $\mathcal{K}_{2}\ni
k\\!-1>\overline{k}\\!-1$, let $i_{k}\in[p]$ be such that
$\lvert(Bx^{k-1})_{i_{k}}\rvert=\lvert Bx^{k-1}\rvert_{\min}$. Since condition
(8) implies that ${\rm supp}(Bx^{k-1})={\rm supp}(B\overline{x}^{k-1})$ for
each $k-1\in\mathcal{K}_{2}$, we have
$\lvert(B\overline{x}^{k-1})_{i_{k}}\rvert\geq\lvert
B\overline{x}^{k-1}\rvert_{\rm min}$. Thus, for each $\mathcal{K}_{2}\ni
k-1>\overline{k}-1$,

$\displaystyle\|Bx^{k-1}-B\overline{x}^{k-1}\|$
$\displaystyle\geq\lvert(Bx^{k-1})_{i_{k}}-(B\overline{x}^{k-1})_{i_{k}}\rvert\geq\lvert(B\overline{x}^{k-1})_{i_{k}}\rvert-\lvert(Bx^{k-1})_{i_{k}}\rvert$
$\displaystyle\geq\lvert B\overline{x}^{k-1}\rvert_{\min}-\lvert
Bx^{k-1}\rvert_{\min}.$

Recall that $\lvert B\overline{x}^{k-1}\rvert_{\min}\geq\nu$ for all
$k\in\mathbb{N}$ by Proposition 4.1 (ii). Together with the last inequality
and $\|Bx^{k-1}-B\overline{x}^{k-1}\|<\frac{\nu}{4}$, for each
$\mathcal{K}_{2}\ni k-1>\overline{k}-1$, we have $\lvert
Bx^{k-1}\rvert_{\min}\geq\frac{3\nu}{4}$. For each $\mathcal{K}_{2}\ni
k-1>\overline{k}-1$, let $j_{k}\in[p]$ be such that
$\lvert(Bx^{k})_{j_{k}}\rvert=\lvert Bx^{k}\rvert_{\min}$. By Remark 4.1 (ii),
${\rm supp}(Bx^{k})\subset{\rm supp}(Bx^{k-1})$ for each
$k-1\in\mathcal{K}_{2}$, which along with $j_{k}\in{\rm supp}(Bx^{k})$ implies
that $\lvert(Bx^{k-1})_{j_{k}}\rvert\geq\lvert Bx^{k-1}\rvert_{\min}$. Thus,
for each $\mathcal{K}_{2}\ni k-1>\overline{k}-1$,

$\displaystyle\|Bd^{k-1}\|$
$\displaystyle\geq\|Bx^{k}-Bx^{k-1}\|\geq\lvert(Bx^{k-1})_{j_{k}}-(Bx^{k})_{j_{k}}\rvert$
$\displaystyle\geq\lvert(Bx^{k-1})_{j_{k}}\rvert-\lvert(Bx^{k})_{j_{k}}\rvert\geq\lvert
Bx^{k-1}\rvert_{\min}-\lvert Bx^{k}\rvert_{\min},$

which together with $\|Bd^{k-1}\|\leq\frac{\nu}{4}$ and $\lvert
Bx^{k-1}\rvert_{\min}\geq\frac{3\nu}{4}$ implies that $\lvert
Bx^{k}\rvert_{\min}\geq\frac{\nu}{2}$.

Claim 2: $T_{k}=\overline{T}_{k}$ for $k>\overline{k}$. From the above
arguments, $\|Bx^{k}\\!-B\overline{x}^{k}\|\leq\frac{\nu}{4}$ for
$k>\overline{k}$. If $i\in T_{k}$, then
$\lvert(B\overline{x}^{k})_{i}\rvert\geq\lvert(Bx^{k})_{i}\rvert-\frac{\nu}{4}\geq\frac{\nu}{4}$,
where the second inequality is using $\lvert
Bx^{k}\rvert_{\min}\\!>\\!\frac{\nu}{2}$ by Claim 1. This means that
$i\in\overline{T}_{k}$, so $T_{k}\subset\overline{T}_{k}$. Conversely, if
$i\in\overline{T}_{k}$, then
$\lvert(Bx^{k})_{i}\rvert\geq\lvert(B\overline{x}^{k})_{i}\rvert-\frac{\nu}{4}\geq\frac{3\nu}{4}$,
so $i\in T_{k}$ and $\overline{T}_{k}\subset T_{k}$. Thus,
$T_{k}=\overline{T}_{k}$ for $k>\overline{k}$.

Claim 3: $T_{k}=T_{k+1}$ for $k>\overline{k}$. If $k\in\mathcal{K}_{1}$, the
result follows directly by the result in Claim 2. If $k\in\mathcal{K}_{2}$,
from the proof of Claim 1,
$\|Bx^{k}\\!-Bx^{k+1}\|\leq\|Bd^{k}\|\leq\frac{\nu}{4}$ for all
$k>\overline{k}$. Then, if $i\in T_{k}$,
$\lvert(Bx^{k+1})_{i}\rvert\geq\lvert(Bx^{k})_{i}\rvert-\frac{\nu}{4}\geq\frac{\nu}{4}$,
where the second inequality is using $\lvert
Bx^{k}\rvert_{\min}\\!>\\!\frac{\nu}{2}$ by Claim 1. This implies that $i\in
T_{k+1}$ and $T_{k}\subset T_{k+1}$. Conversely, if $i\in T_{k+1}$, then
$\lvert(Bx^{k})_{i}\rvert\geq\lvert(Bx^{k+1})_{i}\rvert-\frac{\nu}{4}\geq\frac{\nu}{4}$.
Hence, $i\in T_{k}$ and $T_{k+1}\subset T_{k}$.

From Claim 2 and Claim 3, there exists $T\subset[p]$ such that
$T_{k}=\overline{T}_{k}=T$ for $k>\overline{k}$. Using the similar arguments,
we can also prove that there exists $S\subset[n]$ such that
$S_{k}=\overline{S}_{k}=S$ for all $k>\overline{k}$ (if necessary increasing
$\overline{k}$).

Pick any $x^{*}\in\Gamma(x^{0})$. Let $\\{x^{k}\\}_{k\in K}$ be a subsequence
such that $\lim_{K\ni k\to\infty}x^{k}=x^{*}$. By the above proof, for all
sufficiently large $k\in K$, $\lvert Bx^{k}\rvert_{\min}\geq\frac{\nu}{2}$ and
$\lvert x^{k}\rvert_{\min}\geq\frac{\nu}{2}$, which implies that $\lvert
Bx^{*}\rvert_{\min}\geq\frac{\nu}{2}$ and $\lvert
x^{*}\rvert_{\min}\geq\frac{\nu}{2}$. The results ${\rm supp}(Bx^{*})=T$ and
${\rm supp}(x^{*})=S$ can be obtained by a proof similar to Claim 3. The proof
is completed. $\Box$

By Proposition 5.1, $k\in\mathcal{K}_{2}$ for all $k>\overline{k}$, i.e., the
sequence $\\{x^{k+1}\\}_{k>\overline{k}}$ is generated by the Newton step.
This means that $\\{x^{k+1}\\}_{k>\overline{k}}$ is identical to the one
generated by the inexact projected regularized Newton method starting from
$x^{\overline{k}+1}$. Also, since $\Pi_{k}=\Pi_{*}:=\Pi_{\overline{k}+1}$ for
all $k>\overline{k}$, Algorithm 2 finally reduces to the inexact projected
regularized Newton method for solving

$\min_{x\in\mathbb{R}^{n}}\ \phi(x):=f(x)+\delta_{\Pi_{*}}(x),$ (47)

which is a minimization problem of function $f$ over polyhedron $\Pi_{*}$,
much simpler than the original problem (1). Consequently, the global
convergence and local convergence rate analysis of PGiPN for model (1) boils
down to analyzing those of the inexact projected regularized Newton method for
(47). The rest of this section is devoted to this. Unless otherwise stated,
the notation $\overline{k}$ in the sequel is always the same as that of
Proposition 5.1. In addition, we require the assumption that $\nabla^{2}\\!f$
is locally Lipschitz continuous on $\Gamma(x^{0})$, where $\Gamma(x^{0})$ is
defined in Lemma 5.1 (iii).

###### Assumption 5.1

$\nabla^{2}\\!f$ is locally Lipschitz continuous on an open set
$\mathcal{O}\supset\Gamma(x^{0})$.

Assumption 5.1 is very standard when analyzing the convergence behavior of
Newton-type method. The following lemma reveals that under this assumption,
the step size $\alpha_{k}$ in Newton step takes $1$ when $k$ is sufficiently
large. Since the proof is similar to that of [26, Lemma B.1], the details are
omitted here.

###### Lemma 5.2

Suppose that Assumption 5.1 holds. Then $\alpha_{k}=1$ for sufficiently large
$k$.

Notice that $\Pi_{*}$ is a polyhedron, which can be expressed as

$\Pi_{*}=\big{\\{}x\in\mathbb{R}^{n}\ \lvert\
B_{T_{\overline{k}+1}^{c}\cdot}x=0,\ x_{S_{\overline{k}+1}^{c}}=0,\ x\geq l,\
-x\geq-u\big{\\}}.$ (48)

For any $x\in\mathbb{R}^{n}$, we define multifunction
$\mathcal{A}:\mathbb{R}^{n}\rightrightarrows[2n]$ as

$\mathcal{A}(x):=\\{i\ \lvert\ x_{i}=l_{i}\\}\cup\\{i+n\ \lvert\
x_{i}=u_{i}\\}.$

Clearly, for $x\in\Pi_{*}$, $\mathcal{A}(x)$ is the active set of constraint
$\Pi_{*}$ at $x$. To prove the global convergence for PGiPN, we first show
that $\\{\mathcal{A}(x^{k})\\}_{k\in\mathbb{N}}$ remains stable when $k$ is
sufficiently large, under the following non-degeneracy assumption.

###### Assumption 5.2

For all $x^{*}\in\Gamma(x^{0})$, $0\in\nabla f(x^{*})+{\rm
ri}(\mathcal{N}_{\Pi_{*}}(x^{*})).$

It follows from Proposition 2.1 and Lemma 5.1 (iii) that for each
$x^{*}\in\Gamma(x^{0})$, $x^{*}$ is a stationary point of $F$, which together
with Proposition 5.1 and Lemma 2.1 (i) yields that $0\in\nabla
f(x^{*})+\mathcal{N}_{\Pi_{*}}(x^{*})$, so that Assumption 5.2 substantially
requires that $-\nabla f(x^{*})$ does not belong to the relative
boundary555For $\Xi\subset\mathbb{R}^{n}$, the set difference ${\rm
cl}(\Xi)\backslash{\rm ri}(\Xi)$ is called the relative boundary of $\Xi$, see
[33, p. 44]. of $\mathcal{N}_{\Pi_{*}}(x^{*})$. In the next lemma, we prove
that under Assumptions 5.1-5.2, $\mathcal{A}(x^{k})=\mathcal{A}(x^{k+1})$ for
sufficiently large $k$.

###### Lemma 5.3

Let $\\{x^{k}\\}_{k\in\mathbb{N}}$ be the sequence generated by Algorithm 2.
Suppose that Assumptions 5.1-5.2 hold. Then, there exist
$\mathcal{A}^{*}\subset[2n]$ and a closed and convex cone
$\mathcal{N}^{*}\subset\mathbb{R}^{n}$ such that
$\mathcal{A}(x^{k})=\mathcal{A}^{*}$ and
$\mathcal{N}_{\Pi_{*}}(x^{k})=\mathcal{N}^{*}$ for sufficiently large $k$.

Proof: We complete the proof via the following two claims.

Claim 1: We first claim that $\lim_{k\rightarrow\infty}\|{\rm
proj}_{\mathcal{T}_{\Pi_{*}}(x^{k})}(-\nabla\\!f(x^{k}))\|=0.$ Since $\Pi_{*}$
is polyhedral, for any $x\in\Pi_{*}$, $\mathcal{T}_{\Pi_{*}}(x)$ and
$\mathcal{N}_{\Pi_{*}}(x)$ are closed and convex cones, and
$\mathcal{T}_{\Pi_{*}}(x)$ is polar to $\mathcal{N}_{\Pi_{*}}(x)$, which
implies that when $k$ is sufficiently large, $z={\rm
proj}_{\mathcal{T}_{\Pi_{*}}(x^{k})}(z)+{\rm
proj}_{\mathcal{N}_{\Pi_{*}}(x^{k})}(z)$ holds for any $z\in\mathbb{R}^{n}$.
Then, for all sufficiently large $k$,

$\|{\rm
proj}_{\mathcal{T}_{\Pi_{*}}(x^{k})}(-\nabla\\!f(x^{k}))\|=\|-\\!\nabla\\!f(x^{k})\\!-\\!{\rm
proj}_{\mathcal{N}_{\Pi_{*}}(x^{k})}(-\nabla\\!f(x^{k}))\|={\rm
dist}(0,\partial\phi(x^{k})).$

Therefore, it suffices to prove that $\lim_{k\rightarrow\infty}{\rm
dist}(0,\partial\phi(x^{k}))=0$. By Proposition 5.1, equation (14), Remark 4.1
(iv) and Lemma 5.1 (ii), there exists $\\{\zeta_{k}\\}_{k>\overline{k}}$ with
$\lim_{k\rightarrow\infty}\|\zeta_{k}\|=0$ such that $0\in\nabla
f(x^{k})+G_{k}d^{k}+\zeta_{k}+\mathcal{N}_{\Pi_{k}}(x^{k}+d^{k})$ for each
$k\geq\overline{k}$, which implies that $\nabla f(x^{k}+d^{k})-\nabla
f(x^{k})-G_{k}d^{k}-\zeta_{k}\in\partial\phi(x^{k}+d^{k})$ for each
$k\geq\overline{k}$. This together with Lemma 5.1 (ii) and the continuity of
$\nabla f$ implies that $\lim_{k\rightarrow\infty}{\rm
dist}(0,\partial\phi(x^{k}+d^{k}))=0$. Thus, by Lemma 5.2 we obtain the
desired result.

Claim 2: We now show that $\mathcal{A}(x^{k})\subset\mathcal{A}(x^{k+1})$ for
all sufficiently large $k$. If not, there exists $K\subset\mathbb{N}$ such
that $\mathcal{A}(x^{k})\not\subset\mathcal{A}(x^{k+1})$ for all $k\in K$. If
necessary taking a subsequence, we assume that $\\{x^{k}\\}_{k\in K}$
converges to $x^{*}$. By Lemma 5.1 (ii), $\\{x^{k+1}\\}_{k\in K}$ converges to
$x^{*}$. In addition, from Claim 1 it follows that
$\lim_{k\rightarrow\infty}\|{\rm
proj}_{\mathcal{T}_{\Pi_{*}}(x^{k+1})}(-\nabla\\!f(x^{k+1}))\|=0$. The two
sides along with Assumption 5.2 and [13, Corollary 3.6] yields that
$\mathcal{A}(x^{k+1})=\mathcal{A}(x^{*})$ for all sufficiently large $k\in K$.
Since $\mathcal{A}(x^{k})\subset\mathcal{A}(x^{*})$ for sufficiently large
$k\in K$, we have $\mathcal{A}(x^{k})\subset\mathcal{A}(x^{k+1})$ for
sufficiently large $k\in K$, contradicting to
$\mathcal{A}(x^{k})\not\subset\mathcal{A}(x^{k+1})$ for $k\in K$. The claimed
fact that $\mathcal{A}(x^{k})\subset\mathcal{A}(x^{k+1})$ then follows.

From $\mathcal{A}(x^{k})\subset\mathcal{A}(x^{k+1})$ for all sufficiently
large $k$, $\\{\mathcal{A}(x^{k})\\}_{k\in\mathbb{N}}$ converges to for some
$A^{*}\subset[2n]$ in the sense of Painlev$\acute{e}$-Kuratowski. From the
finiteness of $\mathcal{A}^{*}$, we conclude that
$\mathcal{A}(x^{k})=\mathcal{A}^{*}$ for all sufficiently large $k$. From the
expression of $\Pi_{*}$ in (48) and $\mathcal{A}(x^{k})=\mathcal{A}^{*}$ for
all sufficiently large $k$, we have
$\mathcal{N}_{\Pi_{*}}(x^{k})=\mathcal{N}^{*}$ for all sufficiently large $k$.
$\Box$

Our proof for the global convergence of PGiPN additionally requires the
following assumption.

###### Assumption 5.3

For every sufficiently large $k$, there exists
$\xi_{k}\in\mathcal{N}_{\Pi_{*}}(x^{k})$ such that

$\liminf_{k\rightarrow\infty}\frac{-\langle\nabla\\!f(x^{k})+\xi_{k},d^{k}\rangle}{\|\nabla
f(x^{k})+\xi_{k}\|\|d^{k}\|}>0.$

This assumption essentially requires that for every sufficiently large $k$
there exists one element $\xi_{k}\in\partial\phi(x^{k})$ such that the angle
between $\nabla\\!f(x^{k})+\xi_{k}$ and $d^{k}$ is uniformly larger than
${\pi}/{2}$. For every sufficiently large $k$, since $x^{k}+\alpha
d^{k}\in\Pi_{*}$ for all $\alpha\in[0,1]$, we have
$d^{k}\in\mathcal{T}_{\Pi^{*}}(x^{k})$, which implies that
$\langle\xi^{k},d^{k}\rangle\leq 0$. Together with (39), for every
sufficiently large $k$, the angle between $\nabla\\!f(x^{k})+\xi_{k}$ and
$d^{k}$ is larger than $\pi/2$. This means that it is highly possible for
Assumption 5.3 to hold. Obviously, when $n=1$, it automatically holds.

Next, we show that if $\phi$ is a KL function and Assumptions 5.1-5.3 hold,
the sequence generated by PGiPN is Cauchy and converges to an $L$-stationary
point.

###### Theorem 5.1

Let $\\{x^{k}\\}_{k\in\mathbb{N}}$ be the sequence generated by Algorithm 2.
Suppose that Assumptions 5.1-5.3 hold, and that $\phi$ is a KL function. Then,
$\sum_{k=1}^{\infty}\|x^{k+1}\\!-\\!x^{k}\|<\infty$, and consequently
$\\{x^{k}\\}_{k\in\mathbb{N}}$ converges to an $L$-stationary point of (1).

Proof: If there exists $\widetilde{k}>\overline{k}$ such that
$\phi(x^{\widetilde{k}})=\phi(x^{\widetilde{k}+1})$, then
$F(x^{\widetilde{k}})=F(x^{\widetilde{k}+1})$ by Proposition 5.1, which
together with Lemma 5.1 (i) yields that
$x^{\widetilde{k}}=\overline{x}^{\widetilde{k}}$. Consequently,
$x^{\widetilde{k}}$ meets the termination condition of Algorithm 2, so that
$\\{x^{k}\\}_{k\in\mathbb{N}}$ converges to an $L$-stationary point of (1)
within a finite number of steps. Thus, we only need to consider the case that
$\phi(x^{k})>\phi(x^{k+1})$ for all $k>\overline{k}$. By [12, Lemma 6], there
exist $\varepsilon>0,\eta>0$ and a continuous concave function
$\varphi\in\Upsilon_{\eta}$ such that for all $\overline{x}\in\Gamma(x^{0})$
and $x\in\\{z\in\mathbb{R}^{n}\ \lvert\ {\rm
dist}(z,\Gamma(x^{0}))<\varepsilon\\}\cap[\phi(\overline{x})<\phi<\phi(\overline{x})+\eta]$,
$\varphi^{\prime}(\phi(x)-\phi(\overline{x})){\rm dist}(0,\partial\phi(x))\geq
1$, where $\Upsilon_{\eta}$ is defined in Definition 2.3, and $\Gamma(x^{0})$
is defined in Lemma 5.1 (iii). Clearly, $\lim_{k\rightarrow\infty}{\rm
dist}(x^{k},\Gamma(x^{0}))=0$. Pick any $x^{*}\in\Gamma(x^{0})$. By the
definition of $\phi$, Propositions 4.1 (i) and 5.1, we have
$\lim_{k\to\infty}\phi(x^{k})=\phi(x^{*})$. Then, for $k>\overline{k}$ (if
necessary by increasing $\overline{k}$), $x^{k}\in\\{z\in\mathbb{R}^{n}\
\lvert\ {\rm
dist}(z,\Gamma(x^{0}))<\varepsilon\\}\cap[\phi(x^{*})<\phi<\phi(x^{*})+\eta]$.
Consequently, for all $k>\overline{k}$,

$\varphi^{\prime}(\phi(x^{k})-\phi(x^{*})){\rm
dist}(0,\partial\phi(x^{k}))\geq 1.$ (49)

By Assumption 5.3, there exist $c>0$ and
$\xi_{k}\in\mathcal{N}_{\Pi_{*}}(x^{k})$ such that for all suffciently large
$k$,

$-\langle\nabla f(x^{k})+\xi_{k},d^{k}\rangle>c\|\nabla
f(x^{k})+\xi_{k}\|\|d^{k}\|.$ (50)

From Lemma 5.3 we have
$\mathcal{N}_{\Pi_{*}}(x^{k})=\mathcal{N}_{\Pi_{*}}(x^{k+1})$ for all
$k>\overline{k}$ (by possibly enlarging $\overline{k}$), which implies that
$\xi_{k}\in\mathcal{N}_{\Pi_{*}}(x^{k+1})$. Together with (38), (50) and Lemma
5.2, we have that for all $k>\overline{k}$ (if necessary enlarging
$\overline{k}$),

$\frac{\phi(x^{k})-\phi(x^{k+1})}{{\rm
dist}(0,\partial\phi(x^{k}))}\\!\geq\\!\frac{-\varrho\langle\nabla
f(x^{k})+\xi_{k},d^{k}\rangle}{{\rm
dist}(0,\partial\phi(x^{k}))}\geq\frac{\varrho c\|\nabla
f(x^{k})+\xi_{k}\|\|d^{k}\|}{\|\nabla f(x^{k})+\xi_{k}\|}=\varrho
c\|x^{k+1}-x^{k}\|,$ (51)

where the second inequality follows by $\nabla
f(x^{k})+\xi_{k}\in\partial\phi(x^{k})$ and (50). For each $k$, let
$\Delta_{k}:=\varphi(\phi(x^{k})\\!-\\!\phi(x^{*}))$. From (49), (51) and the
concavity of $\varphi$ on $[0,\eta)$, for all $k>\overline{k}$,

$\displaystyle\Delta_{k}-\Delta_{k+1}$
$\displaystyle=\phi(x^{k})-\phi(x^{k+1})\geq\varphi^{\prime}(\phi(x^{k})\\!-\\!\phi(x^{*}))(\phi(x^{k})\\!-\\!\phi(x^{k+1}))$
$\displaystyle\geq\frac{\phi(x^{k})-\phi(x^{k+1})}{{\rm
dist}(0,\partial\phi(x^{k}))}\geq\varrho c\|x^{k+1}-x^{k}\|.$

Summing this inequality from $\overline{k}$ to any $k>\overline{k}$ and using
$\Delta_{k}\geq 0$ yields that

$\sum_{j=\overline{k}}^{k}\|x^{j+1}\\!-\\!x^{j}\|\leq\frac{1}{\varrho
c}\sum_{j=\overline{k}}^{k}(\Delta_{j}\\!-\\!\Delta_{j+1})=\frac{1}{\varrho
c}(\Delta_{\overline{k}}\\!-\\!\Delta_{k+1})\leq\frac{1}{\varrho
c}\Delta_{\overline{k}}.$

Passing the limit $k\to\infty$ leads to
$\sum_{j=\overline{k}}^{\infty}\|x^{j+1}\\!-\\!x^{j}\|<\infty$. Thus,
$\\{x^{k}\\}_{k\in\mathbb{N}}$ is a Cauchy sequence and converges to $x^{*}$.
It follows from Lemma 5.1 (iii) that $x^{*}$ is an $L$-stationary point of
model (1). The proof is completed. $\Box$

We now focus on the superlinear rate analysis of PGiPN. Denote

$\mathcal{X}^{*}:=\big{\\{}x\in\mathbb{R}^{n}\ \lvert\ 0\in\nabla
f(x)+\mathcal{N}_{\Pi_{*}}(x),\ \nabla^{2}f(x)\succeq 0\big{\\}},$

which we call by the set of second-order stationary points of (47). Based on
this notation, we assume that a local Hölderian error bound condition holds
with $\mathcal{X}^{*}$ in Assumption 5.4. For more introduction on the
Hölderian error bound condition, we refer the interested readers to [30] and
[26].

###### Assumption 5.4

The $q$-subregularity of function $r(x):=x-{\rm proj}_{\Pi_{*}}(x-\nabla
f(x))$ holds at $x^{*}$ for the origin with $\mathcal{X}^{*}$, i.e., there
exist $\varepsilon>0,\ \kappa>0$ and $q\in(0,1]$ such that for all
$x\in\mathbb{B}(x^{*},\varepsilon)\cap\Pi_{*},$ ${\rm dist}(x,r^{-1}(0))={\rm
dist}(x,\mathcal{X}^{*})\leq\kappa\|r(x)\|^{q}.$

Recently, Liu et al. [26] proposed an inexact regularized proximal Newton
method (IRPNM) for solving the problems, consisting of a smooth function and
an extended real-valued convex function, which includes (47) as a special
case. They studied the superlinear convergence rate of IRPNM under Assumptions
5.1 and 5.4. By [37, Lemma 4] and $\overline{\mu}_{k}\in[\mu_{\rm
min},\widetilde{\mu})$, $\|r_{k}(x^{k})\|=O(\|r(x^{k})\|)$ for sufficiently
large $k$. This together with Assumption 5.4 implies that there exists
$\widehat{\kappa}>0$ such that for sufficiently large $k$ with
$x^{k}\in\mathbb{B}(x^{*},\varepsilon)$,

${\rm dist}(x^{k},\mathcal{X}^{*})\leq\widehat{\kappa}\|r_{k}(x^{k})\|^{q}.$
(52)

Recall that PGiPN finally reduces to an inexact projected regularized Newton
method for solving (47). From Lemma 4.1 (iii) and Lemma 5.2, we have

$\Theta_{k}(x^{k+1})-\Theta_{k}(x^{k})\leq 0\ \ {\rm and}\ \
\|R_{k}(x^{k+1})\|\leq\frac{1}{2}\min\\{\|r_{k}(x^{k})\|,\|r_{k}(x^{k})\|^{1+\varsigma}\\},$
(53)

for sufficiently large $k$. Let
$\Lambda_{k}^{i}:=G_{k}^{i}\\!-\\!\nabla^{2}\\!f(x^{k})-b_{1}\|\overline{\mu}_{k}(x^{k}-\overline{x}^{k})\|^{\sigma}I$,
where $G_{k}^{i}$ are those in (34)-(36). Under Assumption 5.4, from [45,
Lemma 4.8] and [26, Lemma 4.4] and the fact that $G_{k}^{1}-G_{k}^{2}\succeq
0$, we have for sufficiently large $k$ with
$x^{k}\in\mathbb{B}(x^{*},\varepsilon)$,

$\max\big{\\{}\lambda_{\min}(\Lambda_{k}^{1}),\lambda_{\min}(\Lambda_{k}^{2}),\lambda_{\min}(\Lambda_{k}^{3})\big{\\}}=O({\rm
dist}(x^{k},\mathcal{X}^{*})),$ (54)

By using (52), (53) and (54), and following a proof similar to [26, Theorem
4.3], we can obtain the following result.

###### Theorem 5.2

Let $\\{x^{k}\\}_{k\in\mathbb{N}}$ be the sequence generated by Algorithm 2.
Suppose that Assumption 5.1 holds, and that $\\{x^{k}\\}_{k\in\mathbb{N}}$
converges to $x^{*}\in\Gamma(x^{0}).$ If Assumption 5.4 holds with
$q\in(\frac{1}{1+\sigma},1]$ at $x^{*}$, then the sequence
$\\{x^{k}\\}_{k\in\mathbb{N}}$ converges to $x^{*}$ with the $Q$-superlinear
convergence rate at order $q(1\\!+\\!\sigma)$.

##  6 Numerical experiments

This section focuses on the numerical experiments of several variants of PGiPN
for solving a fused $\ell_{0}$-norms regularization problem with a box
constraint. We first describe the implementation of Algorithm 2 in Section
6.1. In Section 6.2, we make comparison between model (1) with the least-
squares loss function $f$ and the fused Lasso model (5) by using PGiPN to
solver the former and SSNAL [23] to solve the latter, to highlight the
advantages and disadvantages of our proposed fused $\ell_{0}$-norms
regularization. Among others, the code of SSNAL is available at
https://github.com/MatOpt/SuiteLasso). Finally, in Section 6.3, we present
some numerical results toward the comparison among several variants of PGiPN
and ZeroFPR and PG method for (1) in terms of efficiency and the quality of
the output. The MATLAB code of PGiPN is available at
https://github.com/yuqiawu/PGiPN.

###  6.1 Implementation of Algorithm 2

####  6.1.1 Dimension reduction of (11)

Suppose that $\emptyset\neq S_{k}^{c}:=[n]\backslash S_{k}$. Based on the fact
that every $x\in\Pi_{k}$ satisfies $x_{S_{k}^{c}}=0$, we can obtain an
approximate solution to (11) by solving a problem in a lower dimension.
Specifically, for each $k\in\mathcal{K}_{2}$, write

$H_{k}\\!:=\\!(G_{k})_{S_{k}S_{k}},\,v^{k}\\!:=\\!x^{k}_{S_{k}},\,\nabla\\!f_{S_{k}}(v^{k})\\!=\\![\nabla\\!f(x^{k})]_{S_{k}},\,\widehat{\Pi}_{k}\\!:=\\{v\in\mathbb{R}^{\lvert
S_{k}\rvert}\ \lvert\
\widetilde{B}_{k}v=0,\,l_{S_{k}}\\!\leq\\!v\\!\leq\\!u_{S_{k}}\\},$

where $\widetilde{B}_{k}$ is the matrix obtained by removing the rows of
$B_{T_{k}^{c}S_{k}}$ whose elements are all zero. We turn to consider the
following strongly convex optimization problem,

$\widehat{v}^{k}\approx\mathop{\arg\min}_{v\in\mathbb{R}^{\lvert
S_{k}\rvert}}\Big{\\{}\theta_{k}(v):=f(I_{\cdot
S_{k}}v^{k})+\langle\nabla\\!f_{S_{k}}(v^{k}),v-\\!v^{k}\rangle+\frac{1}{2}(v\\!-v^{k})^{\top}H_{k}(v\\!-\\!v^{k})+\delta_{\widehat{\Pi}_{k}}(v)\Big{\\}}.$
(55)

The following lemma gives a way to find $y^{k}$ satisfying (13)-(14) by
inexactly solving problem (55), whose dimension is much smaller than that of
(11) if $\lvert S_{k}\rvert\ll n$.

###### Lemma 6.1

Let $y^{k}_{S_{k}}=\widehat{v}^{k}$ and $y^{k}_{S_{k}^{c}}=0$. Then,
$\Theta_{k}(y^{k})=\theta_{k}(\widehat{v}^{k})$ and ${\rm
dist}(0,\partial\Theta_{k}(y^{k}))={\rm
dist}(0,\partial\theta_{k}(\widehat{v}^{k}))$. Consequently, the vector
$\widehat{v}^{k}$ satisfies

$\theta_{k}(\widehat{v}^{k})-\theta_{k}(v^{k})\leq 0,\ {\rm
dist}(0,\partial\theta_{k}(\widehat{v}^{k}))\leq\frac{\min\\{\overline{\mu}_{k}^{-1},1\\}}{2}\min\left\\{\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|,\|\overline{\mu}_{k}(x^{k}\\!-\\!\overline{x}^{k})\|^{1+\varsigma}\right\\},$

if and only if the vector $y^{k}$ satisfies the inexact conditions in
(13)-(14).

Proof: The first part is straightforward. We consider the second part. By the
definition of $\Theta_{k}$, ${\rm dist}(0,\partial\Theta_{k}(y^{k}))={\rm
dist}(0,\nabla f(x^{k})+G_{k}(y^{k}-x^{k})+\mathcal{N}_{\Pi_{k}}(y^{k}))$.
Recall that $\Pi_{k}=\\{x\in\Omega\ \lvert\
B_{T_{k}^{c}\cdot}x=0,x_{S_{k}^{c}}=0\\}$. Then,
$\mathcal{N}_{\Pi_{k}}(y^{k})={\rm Range}(B_{T_{k}^{c}\cdot}^{\top})+{\rm
Range}(I^{\top}_{S_{k}^{c}\cdot})+\mathcal{N}_{\Omega}(y^{k})$, and

$\displaystyle{\rm dist}(0,\partial\Theta_{k}(y^{k}))$ $\displaystyle={\rm
dist}\big{(}0,\nabla f(x^{k})+G_{k}(y^{k}-x^{k})+{\rm
Range}(B_{T_{k}^{c}\cdot}^{\top})+{\rm
Range}(I^{\top}_{S_{k}^{c}\cdot})+\mathcal{N}_{\Omega}(y^{k})\big{)}$
$\displaystyle={\rm dist}\big{(}0,\nabla
f_{S_{k}}(v^{k})+H_{k}(\widehat{v}^{k}-v^{k})+{\rm
Range}(B_{T_{k}^{c}S_{k}}^{\top})+\mathcal{N}_{[l_{S_{k}},u_{S_{k}}]}(\widehat{v}^{k})\big{)}$
$\displaystyle={\rm dist}(0,\nabla
f_{S_{k}}(v^{k})+H_{k}(\widehat{v}^{k}-v^{k})+\mathcal{N}_{\widehat{\Pi}_{k}}(\widehat{v}^{k}))={\rm
dist}(0,\theta_{k}(\widehat{v}^{k})),$

where the second equality is using ${\rm
Range}(I^{\top}_{S_{k}^{c}\cdot})=\\{z\in\mathbb{R}^{n}\ |\ z_{S_{k}}=0$}.
$\Box$

####  6.1.2 Acceleration of Algorithm 1

The switch condition in (8) is in general difficult to be satisfied when
$\|Bx^{k}\|_{0}$ or $\|x^{k}\|_{0}$ is large. Consequently, PGiPN is
continuously executing PG steps. This phenomenon is evident in the numerical
experiment of the restoration of blurred image, see Section 6.3.2. To further
accelerate the iterates of Algorithm 2 into the Newton step, we introduce the
following relaxed switch condition:

$\|\lvert{\rm sign}(Bx^{k})\rvert-\lvert{\rm
sign}(B\overline{x}^{k})\rvert\|_{1}\leq\frac{\eta_{1}n}{k}\ \ {\rm and}\ \
\|\lvert{\rm sign}(x^{k})\rvert-\lvert{\rm
sign}(\overline{x}^{k})\rvert\|_{1}\leq\frac{\eta_{2}n}{k},$ (56)

where $\eta_{1},\eta_{2}$ are two nonnegative constants. Following the
arguments similar to those in Lemma 4.1, we have that Algorithm 2 equipped
with (56) is also well defined. Obviously, when $\frac{\eta_{i}n}{k}\geq 1$,
condition (56) allows the supports of $Bx^{k}$ and $B\overline{x}^{k}$ and
$x^{k}$ and $\overline{x}^{k}$ have some difference; when
$\frac{\eta_{i}n}{k}<1(i=1,2)$, condition (56) is identical to (8). This means
that as $k$ grows, Algorithm 2 with relaxed switch condition (56) will finally
reduce to the one with (8). Since our convergence analysis does not specify
the initial point, the asymptotic convergence results also hold for Algorithm
2 with condition (56).

####  6.1.3 Choice of parameters in Algorithm 2

We will test the performance of PGiPN with $G_{k}$ given by $G_{k}^{2}$ in
(35), and PGiPN(r), which is PGiPN with relaxed switch condition (56). We use
Gurobi to solve subproblem (11) with such $G_{k}$, with inexact conditions
(13), (14) controlled by options params.Cutoff and params.OptimalityTol,
respectively. Also, we test PGilbfgs, which is the same as PGiPN, except using
limited-memory BFGS (lbfgs) to construct $G_{k}$. In particular, we form
$G_{k}=B_{k}+b_{1}\|\overline{\mu}_{k}(x^{k}-\overline{x}^{k})\|^{\sigma}$,
with $B_{k}$ given by lbfgs. For solving (11) with such $G_{k}$, we use the
method introduced in [19]. We set the parameters of all the variants of PGiPN
by $\alpha=10^{-8},\ \sigma=\frac{1}{2},\ \varrho=10^{-4},\
\beta=\frac{1}{2},\ \varsigma=\frac{2}{3}.$ We set $b_{1}=10^{-3}$ for PGiPN
and PGiPN(r), and $b_{1}=10^{-8}$ for PGilbfgs.

We compare the numerical performance of our algorithms with those of ZeroFPR
[40] and the PG method [43]. In particular, ZeroFPR uses the quasi-Newton
method to minimize the forward-backward envelope of the objective. The code
package of ZeroFPR is downloaded from http://github.com/kul-forbes/ForBES. We
set “lbfgs” as the solver of ZeroFPR. On the other hand, the iterate steps of
PG are the same as those of PGiPN without the Newton steps, so that we can
check the effect of the additional second-order step on PGiPN. For this
reason, the parameters of PG are chosen to be the same as those involved in PG
Step of PGiPN. We also observe that the sparsity of the output is very
sensitive to $\mu_{k}$ in Algorithm 2. To be fair, as the default setting in
ZeroFPR, in all variants of PGiPN and PG, we set $\mu_{k}=0.95^{-1}L_{1}$ for
all $k\in\mathbb{N}$, where $L_{1}$ is an estimation of the Lipschitz constant
of $\nabla\\!f$ obtained by computing $\|A\|_{2}$ from the following MATLAB
sentences:

opt.issym = 1; opt.tol = 0.001; ATAmap = @(x) A’*A*x; L =
eigs(ATAmap,n,1,‘LM’,opt).

For each solver, we set $x^{0}=0$ and terminate at the iterate $x^{k}$
whenever $k\geq 5000$ or $\overline{\mu}_{k}\|x^{k}-{\rm
prox}_{\overline{\mu}_{k}^{-1}g}(x^{k}-\overline{\mu}_{k}^{-1}\nabla\\!f(x^{k}))\|_{\infty}<10^{-4}$.
All the numerical tests in this section are conducted on a desktop running on
64-bit Windows System with an Intel(R) Core(TM) i7-10700 CPU 2.90GHz and 32.0
GB RAM.

###  6.2 Model comparison with the fused Lasso

This subsection is devoted to the numerical comparison between the fused
$\ell_{0}$-norms regularization problem with a box constraint (FZNS), i.e.,
model (1) with $f=\frac{1}{2}\|A\cdot-b\|^{2}$ and $B=\widehat{B}$ and the
fused Lasso (5). We apply PGiPN to solve FZNS and SSNAL to solve (5). Since
the solved models are different, we only compare the quality of solutions
returned by these two solvers, and will not compare their running time.

Our first empirical study focuses on the ability of regression. For this
purpose, we use a commonly used dataset, prostate data, which can be
downloaded from https://hastie.su.domains/ElemStatLearn/. There are 97
observations and 9 features included in this dataset. This data was used in
[18] to check the performance of square root fused Lasso.

We randomly select 50 observations to form the training set, which composes
$A\in\mathbb{R}^{50\times 8}$. The corresponding responses are represented by
$b\in\mathbb{R}^{50}$. The reminders are left as testing set, which forms
$(\bar{A},\bar{b})$ with $\bar{A}\in\mathbb{R}^{47\times 8}$ and
$\bar{b}\in\mathbb{R}^{47}$. We employ PGiPN to solve FZNS, and SSNAL [23] to
solve the fused Lasso (5), with $(A,b)$ given above, and $l=-1000\times{\bf
1}$, $u=1000\times{\bf 1}$. For each solver, we select $10$ groups of
$(\lambda_{1},\lambda_{2})\in[0.003,400]\times[0.0003,40]$, ensuring that the
outputs exhibit different sparsity levels. We record the sparsity and the
testing error, where the later one is defined as
$\|\bar{A}x^{*}\\!-\\!\bar{b}\|$ with $x^{*}$ being the output. The above
procedure is repeated for 100 randomly constructed $(A,b)$, resulting in a
total of 1000 recorded outputs for each model. All the sparsity pairs
$(\|\widehat{B}x^{*}\|_{0},\|x^{*}\|_{0})$ from PGiPN and SSNAL are recorded
in lines 1, 3 and 5 in Table 1. For each sparsity pair, the mean testing
errors of $\|\bar{A}x^{*}\\!-\\!\bar{b}\|$ for PGiPN and SSNAL corresponding
to the given pair is recorded in lines 2, 4 and 6 in Table 1. Among others,
since the fused Lasso may produce solutions with components being very small
but not equal to $0$, we define $\|y\|_{0}:=\min\\{k\ |\
\sum_{i=1}^{k}|\hat{y}|\geq 0.999\|y\|_{1}\\}$ as in [23], where $\hat{y}$ is
obtained by sorting $y$ in a nonincreasing order, for the outputs of the fused
Lasso. As shown in Table 1, it is evident that when
$(\|\widehat{B}x^{*}\|_{0},\|x^{*}\|_{0})=(6,6)$, the mean testing error for
FZNS is the smallest among all the testing examples. Furthermore, in the
presented 21 comparative experiments, the fused Lasso outperforms FZNS for
only 8 cases. Among these 8 experiments, in 7 cases,
$\|\widehat{B}x^{*}\|_{0}\geq 4$ and $\|x^{*}\|_{0}\geq 6$. This indicates
that our model performs better when the solution is relatively sparse.

Table 1: Mean testing error (FZNS$|$Fused Lasso) of the outputs. ($\|\widehat{B}x^{*}\|_{0},\|x^{*}\|_{0}$) | (2,1) | (3,2) | (3,3) | (3,4) | (3,5) | (3,6) | (3,8)
---|---|---|---|---|---|---|---
Mean testing error |  8.35$|$8.54 |  7.34$|$7.36 |  5.45$|$5.15 |  5.15$|$5.74 |  5.21$|$6.32 |  5.08$|$5.27 |  5.11$|$5.70
($\|\widehat{B}x^{*}\|_{0},\|x^{*}\|_{0}$) | (4,4) | (4,5) | (4,6) | (4,7) | (4,8) | (5,5) | (5,6)
Mean testing error |  5.07$|$5.52 |  5.31$|$5.86 |  5.49$|$4.99 |  5.25$|$4.97 |  5.33$|$4.78 |  5.10$|$5.48 |  5.74$|$5.38
($\|\widehat{B}x^{*}\|_{0},\|x^{*}\|_{0}$) | (5,7) | (5,8) | (6,6) | (6,7) | (6,8) | (7,7) | (7,8)
Mean testing error |  5.46$|$5.58 |  5.35$|$ 5.19 |  4.41$|$ 5.26 |  5.34$|$4.95 |  5.25$|$5.34 |  5.03$|$5.22 |  5.24$|$5.22

Our second numerical study is to evaluate the classification ability of the
two models using the TIMIT database. As introduced in Section 1, the TIMIT
database is a widely used resource for research in speech recognition.
Following the approach described in [21], we compute a log-periodogram from
each speech frame, which is one of the several widely used methods to generate
speech data in a form suitable for speech recognition. Consequently, the
dataset comprises 4509 log-periodograms of length 256 (frequency). It was
highlighted in [21] that distinguishing between “aa” and “ao” is particularly
challenging. Our aim is to classify these sounds using FZNS and the fused
Lasso with $\lambda_{2}=0,\ l=-{\bf 1}$ and $u={\bf 1}$, or in other words,
the zero order variable fusion (3) plus a box constraint and the first order
variable fusion (4).

In TIMIT, the numbers of phonemes labeled “aa” and “ao” are 695 and 1022,
respectively. As in [21], we use the first 150 frequencies of the log-
periodograms because the remaining 106 frequencies do not appear to contain
any information. We randomly select $m_{1}$ samples labeled “aa” and $m_{2}$
samples labeled “ao” as training set, which together with their labels form
$A\in\mathbb{R}^{m\times n}$ and $b\in\mathbb{R}^{m}$, with $m=m_{1}+m_{2},\
n=150$, where $b_{i}=1$ if $A_{i\cdot}$ is labeled as “aa”, and $b_{i}=2$
otherwise. The rest of dataset is left as the testing set, which forms
$\bar{A}\in\mathbb{R}^{(1717-m)\times n},\bar{b}^{1717-m}$, with
$\bar{b}_{i}=1$ if $\bar{A}_{i\cdot}$ is labeled as “aa” and $\bar{b}_{i}=2$
otherwise. For $(A,b)$, given 10 $\lambda_{1}$’s randomly selected within
$[2\times 10^{-5},300]$ such that the sparsity of the outputs
$\|\widehat{B}x^{*}\|_{0}$ spans a wide range. If $\bar{A}_{i\cdot}x^{*}\leq
1.5$, this phoneme is classified as “aa” and hence we set $\hat{b}_{i}=1$;
otherwise, $\hat{b}_{i}=2$. If $\hat{b}_{i}\neq\bar{b}_{i}$, $A_{i\cdot}$ is
regarded as failure in classification. Then the error rate of classification
is given by $\frac{\|\bar{b}-\hat{b}\|_{1}}{1717-m}$. We record both
$\|\widehat{B}x^{*}\|_{0}$ and the error rate of classification.

The above procedure is repeated for $30$ groups of randomly generated $(A,b)$,
resulting in 300 outputs for each solver. The four figures in Figure 1 present
$\|\widehat{B}x^{*}\|_{0}$ and the error rate for each output, with 4
different choices of $(m_{1},m_{2})$. We can see that, for each figure the
output with the smallest error rate is always achieved by the fused
$\ell_{0}$-norms regularization model. It is apparent that in general, FZNS
performs better than the fused Lasso when $\|\widehat{B}x^{*}\|_{0}\leq 30$,
while the mean error rate of the fused Lasso is lower than that of FZNS when
$\|\widehat{B}x^{*}\|_{0}\geq 60$. This phenomenon is especially evident when
$m_{1}$ and $m_{2}$ are small.

(a) $m_{1}=25,m_{2}=50$

(b) $m_{1}=50,m_{2}=100$

(c) $m_{1}=100,m_{2}=200$

(d) $m_{1}=200,m_{2}=400$

Figure 1: $\|\widehat{B}x^{*}\|_{0}$ and the classification error rate for the
outputs from FZNS and the fused Lasso under different $m_{1},m_{2}$.

Based on the results of these two empirical studies, we deduce that the fused
$\ell_{0}$-norms regularization tends to outperform the fused Lasso
regularization model when the output is sufficiently sparse. However, it is
important to note that the numerical performance of the fused $\ell_{0}$-norms
regularization is not stable if the output is not sparse, especially when the
number of observations is small, which suggests that when employing the fused
$\ell_{0}$-norms regularization, careful consideration should be given to
selecting an appropriate penalty parameter. Moreover, due to the fact that for
some optimal solution $x^{*}$ of the fused Lasso regularization problem,
$|\widehat{B}x^{*}|_{\min}$ and $|x^{*}|_{\min}$ may be very small but not
equal to zero, which leads to a difficulty in interpreting what the outputs
mean in the real world application. This also well matches the statements in
[21] that the $\ell_{0}$-norm variable fusion produces simpler estimated
coefficient vectors.

###  6.3 Comparison with ZeroFPR and PG

This subsection focuses on the comparison among several variants of PGiPN,
ZeroFPR and PG, in terms of efficiency and the quality of the outputs.

####  6.3.1 Classification of TIMIT

The experimental data used in this part is the TIMIT dataset, the one in
Section 6.2. To test the performance of the algorithms on (1) with nonconvex
$f$, we consider solving model (1) with $f=\sum_{i=1}^{m}\
\log\Big{(}1+\frac{(A\cdot-b)_{i}}{\nu}\Big{)}$, $B=\widehat{B}$, $l=-{\bf 1}$
and $u={\bf 1}$, where $A\in\mathbb{R}^{m\times n}$ represents the training
data and $b\in\mathbb{R}^{m}$ is the vector of the corresponding labels. It is
worth noting that the loss function is nonconvex, and it was claimed in [2]
that this loss function is effective to process data denoised by heavy-tailed
Student’s $t$-noise.

Following the approach in Section 6.2, we use the first 150 frequencies of the
log-periodograms. For the training set, we arbitrarily selecte 200 samples
labeled as “aa” and 400 samples labeled as “ao”. These samples, along with
their corresponding labels, form the matrices $A\in\mathbb{R}^{m\times n}$ and
$b\in\mathbb{R}^{m}$, with dimensions $m=600$ and $n=150$. The remaining
samples are designated as the testing set. Given a series of nonnegative
$\lambda_{c}$, we set $\lambda_{1}=\lambda_{c}\times
10^{-7}\|A^{\top}b\|_{\infty}$ and $\lambda_{2}=0.1\lambda_{1}$. We employ
four solvers: PGiPN, PGilbfgs, PG, and ZeroFPR. Subsequently, we record the
CPU time and the error rate of classification on the testing set. This
experimental procedure is repeated for a total of 30 groups of $(A,b)$, and
the mean CPU time and error rate are recorded for each $\lambda_{c}$,
presented in Figure 2. Motivated by the experiment in Section 6.2, we also
draw Figure 3, recording $\|\widehat{B}x^{*}\|_{0}$ and the error rate for all
the tested cases for four solvers.

(a) $\lambda_{c}$-$\log({\rm time(seconds)})$ plot

(b) $\lambda_{c}$-error rate plot

Figure 2: Mean of the cpu time and the error rate on 30 examples for four
solvers Figure 3: Scatter figure for all tested examples, recording the
relationship of sparsity ($\|\widehat{B}x^{*}\|_{0}$) and the error rate of
classification.

We see from Figure 2(a) that in terms of efficiency, PGiPN is always the best
one, more than ten times faster than the other three solvers. The reason is
that the other three solvers depend heavily on the proximal mapping of $g$,
and its computation is a little time-consuming, which reflects the advantage
of the projected regularized Newton steps in PGiPN. From Figure 2(b), when
$\lambda_{c}=1$, PGiPN reaches the smallest mean error rate among four solvers
for 10 $\lambda_{c}$’s. When $\lambda_{c}$ is large ($>0.4$), PGiPN and
PGilbfgs tend to outperform ZeroFPR and PG. Moreover, when $\lambda_{c}$ is
small $(<0.1)$, the solutions returned by PG have the best error rate among
four solvers. This is because $\widehat{B}x^{*}$ produced by PG is sparser
than those of the other three solvers under the same $\lambda_{c}$, which can
be observed from Figure 3. For small $\lambda_{c}$, the solutions by the other
three solvers are not sparse, leading to high error rate of classification.

####  6.3.2 Recovery of blurred images

Let $\overline{x}\in\mathbb{R}^{n}$ with $n=256^{2}$ be a vector obtained by
vectorizing a $256\\!\times\\!256$ image “cameraman.tif” in MATLAB, and be
scaled such that all the entries belong to $[0,1]$. Let
$A\in\mathbb{R}^{n\times n}$ be a matrix representing a Gaussian blur operator
with standard deviation $4$ and a filter size of $9$, and the vector
$b\in\mathbb{R}^{m}$ represents a blurred image obtained by adding Gauss noise
$e\sim\mathcal{N}(0,\epsilon)$ with $\epsilon>0$ to $A\overline{x}$, i.e.,
$b=A\overline{x}+e$. We apply model (1) with $f=\frac{1}{2}\|A\cdot-b\|^{2}$,
$B=\widehat{B}$, $l=0$ and $u={\bf 1}$, to restore the blurred images. We test
five solvers, which are PGiPN, PGiPN(r), PGilbfgs, ZeroFPR and PG. For
PGiPN(r), we set $\eta_{1}=0.01,\eta_{2}=0.01$ in (56). For all these five
solvers, we employ
$\lambda_{1}=\lambda_{2}=0.0005\times\|A^{\top}b\|_{\infty}$. Under different
$\epsilon$’s, we compare the performance of these five solvers in terms of
required iterations (Iter), cpu time (Time), $F(x^{*})$ (Fval),
$\|x^{*}\|_{0}$ (xNnz), $\|\widehat{B}x^{*}\|_{0}$ (BxNnz) and the highest
peak signal-to-noise ratio (PSNR), where ${\rm
PSNR}:=10\log_{10}\left(\frac{n}{\|\overline{x}-x^{*}\|^{2}}\right).$ In
particular, to check the effect of the Newton step, we record the iterations
(or time) in the form $M(N)$, where $M$ means the total iterations (or time)
and $N$ means the iterations (or time) in Newton step. PSNR measures the
quality of the restored images. The higher PSNR, the better the quality of
restoration. Table 2 presents the numerical results.

Table 2: Numerical comparison of six solvers on recovery of blurred image with $\lambda_{1}=\lambda_{2}=0.0005\|A^{\top}b\|_{\infty}$ Noise |  |  | PGiPN | PGiPN(r) | PGilbfgs | PG | ZeroFPR
---|---|---|---|---|---|---|---
$\epsilon=0.01$ |  | Iter | 379(3) | 123(6) | 529(31) | 796 | 361
| Time | 1.70e3(10.2) | 5.61e2(21.4) | 2.39e3(6.0) | 3.43e3 | 2.35e4
| Fval | 37.88 | 37.95 | 37.88 | 37.88 | 37.77
| xNnz | 63858 | 63805 | 63858 | 63858 | 63717
| BxNnz | 5767 | 5995 | 5776 | 5779 | 5834
| psnr | 25.90 | 25.77 | 25.90 | 25.90 | 25.91
$\epsilon=0.02$ |  | Iter | 281(4) | 109(6) | 457(38) | 853 | 286
| Time | 1.20e3(8.8) | 4.74e2(13.2) | 1.95e3(7.4) | 3.62e3 | 1.82e4
| Fval | 45.98 | 46.05 | 45.98 | 45.98 | 45.83
| xNnz | 63495 | 63440 | 63495 | 63495 | 63350
| BxNnz | 6098 | 6320 | 6098 | 6099 | 6143
| psnr | 25.41 | 25.23 | 25.42 | 25.42 | 25.33
$\epsilon=0.03$ |  | Iter | 234(3) | 94(3) | 325(15) | 717 | 332
| Time | 9.8e2(6.5) | 3.98e2(6.5) | 1.36e3(6.7) | 2.97e3 | 1.88e4
| Fval | 60.26 | 60.34 | 60.26 | 60.26 | 60.02
| xNnz | 63006 | 62944 | 63006 | 63006 | 62800
| BxNnz | 6594 | 6844 | 6597 | 6592 | 6710
| psnr | 24.90 | 24.67 | 24.90 | 24.90 | 24.76
$\epsilon=0.04$ |  | Iter | 255(3) | 78(5) | 360(19) | 526 | 230
| Time | 1.04e3(6.2) | 3.37e2(20.0) | 1.45e3(3.9) | 2.18e3 | 1.11e4
| Fval | 77.82 | 77.87 | 77.82 | 77.82 | 77.44
| xNnz | 62103 | 62002 | 62104 | 62104 | 61853
| BxNnz | 7267 | 7553 | 7268 | 7271 | 7427
| psnr | 24.20 | 23.85 | 24.20 | 24.20 | 24.00
$\epsilon=0.05$ |  | Iter | 263(3) | 76(11) | 389(29) | 688 | 168
| Time | 1.05e3(10.3) | 3.30e2(28.5) | 1.55e3(6.3) | 2.71e3 | 5.91e3
| Fval | 99.65 | 99.72 | 99.65 | 99.65 | 98.93
| xNnz | 61376 | 61286 | 61381 | 61381 | 60963
| BxNnz | 7955 | 8283 | 7955 | 7956 | 8240
| psnr | 23.36 | 23.00 | 23.37 | 23.37 | 22.87

From Table 2, PGiPN(r) always performs the best in terms of efficiency, which
verifies the effectiveness of the acceleration scheme proposed in Section
6.1.2. PGiPN is faster than PGilbfgs, and PGilbfgs is faster than PG,
supporting the effective acceleration of the Newton steps. However, ZeroFPR is
the most time-consuming, even worse than PG, a pure first-order method. The
reason is that ZeroFPR requires more line searches, and each line search
involves a computation of the proximal mapping of $g$, which is expensive (2-5
seconds).

Despite the superiority of efficiency, the solutions yielded by PGiPN(r) is
not good. We also observe that $\|\widehat{B}x^{*}\|_{0}$ of PGiPN(r) is a
little higher than those of PGiPN, PGilbfgs and PG, because PGiPN(r) runs few
PG steps, so that its structured sparsity is not well reduced. Moreover, the
PSNR is closely related to $\|\widehat{B}x^{*}\|_{0}$, and this leads to the
weakest performance of PGiPN(r) in terms of PSNR. On the other hand, although
ZeroFPR always outputs solutions with the smallest objective value, its PSNR
is not as good as the objective value. The performance of PGiPN, PGilbfgs and
PG in terms of the objective value and PSNR are quite similar. Taking the
efficiency and the quality of the output into consideration, we conclude that
PGiPN is the best solver for this test.

##  7 Conclusions