2021

[5]Jein-Shan Chen

1]Department of Mathematics, National Cheng Kung University, Tainan, 701,
Taiwan

2]Institute of Statistical Science, Academia Sinica, Taipei, 11677, Taiwan

3]Center for Advanced Intelligence Project, RIKEN, Tokyo, 103-0027, Japan

4]Department of Creative Informatics, Graduate School of Information Science
and Technology, The University of Tokyo, Tokyo, 113-8656, Japan; Center for
Advanced Intelligence Project, RIKEN, Tokyo, 103-0027, Japan

[5]Department of Mathematics, National Taiwan Normal University, Taipei,
11677, Taiwan

# Unified Smoothing Approach for Best Hyperparameter Selection Problem Using a
Bilevel Optimization Strategy

Chieu Thanh Nguyen<EMAIL_ADDRESS>Jan Harold Alcantara
<EMAIL_ADDRESS>Takayuki Okuno<EMAIL_ADDRESS>Akiko
Takeda<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  *

###### Abstract

Strongly motivated from use in various fields including machine learning, the
methodology of sparse optimization has been developed intensively so far.
Especially, the recent advance of algorithms for solving problems with
nonsmooth regularizers is remarkable. However, those algorithms suppose that
weight parameters of regularizers, called hyperparameters hereafter, are pre-
fixed, and it is a crucial matter how the best hyperparameter should be
selected. In this paper, we focus on the hyperparameter selection of
regularizers related to the $\ell_{p}$ function with $0<p\leq 1$ and consider
to apply the bilevel programming strategy, wherein we need to solve a bilevel
problem, whose lower-level problem is a nonsmooth, possibly nonconvex and non-
Lipschitz. Recently, for solving a bilevel problem for hyperparameter
selection of the pure $\ell_{p}\ (0<p\leq 1)$ regularizer Okuno et al has
discovered the new necessary optimality conditions, called SB(scaled
bilevel)-KKT conditions, and further proposed a smoothing-type algorithm using
a certain specific smoothing function.

Our contribution is to propose a unified smoothing approach using smoothing
functions that belong to the Chen-Mangasarian class, and then prove that
generated iteration points accumulate at SB-KKT points under milder
assumptions than Okuno et al. Another contribution is that our approach and
analysis are applicable to a wider class of regularizers than Okuno et al.
Numerical comparisons demonstrate which smoothing functions work well for
hyperparameter optimization via bilevel optimization approach.

###### keywords:

hyperparameter learning, smoothing functions, bilevel optimization

##  1 Introduction

A learning algorithm in machine learning usually involves solving the
unconstrained optimization problem

$\min_{\omega\in\Re^{n}}g(\omega)+\sum_{i=1}^{r}\lambda_{i}R_{i}(\omega),$ (1)

where $\lambda=(\lambda_{1},\dots,\lambda_{r})$ is called a hyperparameter,
whose value is decided prior to implementation of the learning algorithm.
Here, $R_{i},g:\Re^{n}\to\Re$, $i=2,\dots,r$ are twice continuously
differentiable functions, and
$R_{1}(\omega):=\sum_{i=1}^{n}\psi(\lvert\omega_{i}\rvert^{p})$ $(0<p\leq 1)$
with $\psi$ satisfying the following assumption:

  * (A)

$\psi:[0,\infty)\to\Re$ is twice continuously differentiable on $[0,\infty)$
and there exist two positive constants $\alpha,\beta$ such that $\forall
t\in[0,\infty)$

$0<\psi^{\prime}(t)\leq\alpha~{}~{}\mbox{and}~{}~{}-\beta\leq\psi^{\prime\prime}(t)\leq
0.$

In this manuscript, we make (A) our blanket assumption on $\psi$. It is well-
known that the function $R_{1}(\omega)$ is nonsmooth, nonconvex, even non-
Lipschitz when $p\in(0,1)$. There are many penalty functions often used in
statistics and signal reconstruction satisfying Assumption (A) (see Appendix
A).

For notation purposes, we denote

$G(\omega,\bar{\lambda}):=g(\omega)+\sum_{i=2}^{r}\lambda_{i}R_{i}(\omega),$

with $\bar{\lambda}:=(\lambda_{2},\dots,\lambda_{r})^{T}\in\Re^{r-1}$. Then
problem (1) can be rewritten as

$\min_{\omega\in\Re^{n}}G(\omega,\bar{\lambda})+\lambda_{1}R_{1}(\omega).$ (2)

The problem of finding optimal values of the hyperparameters can be
accomplished using grid search and Bayesian optimization BB12 ; WHZMF16 . This
paper, on the other hand, is devoted to a bilevel optimization strategy to
find the best hyperparameter.

In particular, we focus on the bilevel nonsmooth programming problem

$\begin{array}[]{cc}\displaystyle{\min_{\omega^{*}_{\lambda},\lambda}}&f(\omega^{*}_{\lambda})\\\
{\rm
s.t}&\omega^{*}_{\lambda}\in\underset{\omega\in\Re^{n}}{\mbox{argmin}}\;G(\omega,\bar{\lambda})+\lambda_{1}R_{1}(\omega)\\\
&\lambda\geq 0,\end{array}$ (3)

where $f:\Re^{n}\to\Re$ is continuously differentiable. Problem (2) which
appears in the constraint set of (3) is called the lower-level problem, and
the minimization of $f$ is called the upper-level problem.

Bilevel optimization problems were introduced by Bracken and McGill BM73 . The
reader is referred to D03 ; CMS07 ; SMD18  for a survey of methods for solving
the bilevel optimization problem and their applications. Efforts have been put
forth by many researchers in the past few decades to use bilevel optimization
strategy to the problem of finding the best hyperparameter values. In
particular, BHJKP06 ; BKJP08  focused on a bilevel support-vector regression
(SVR) problem where the lower-level optimization problem is cast as a convex
quadratic program. MBB09 ; MBB11  proposed a bilevel cross-validation program
for support-vector machine (SVM), where the upper-level problem is convex and
nonsmooth, while the lower-level problem is differentiable. ORBP16  used
gradient-based methods for the bilevel optimization problem with nonsmooth
convex lower-level problem (for example, sparse models based on the
$\ell_{1}$-norm). However, BHJKP06 ; MBB09 ; ORBP16  only proposed the
algorithms to solve the bilevel optimization without providing any theoretical
analysis. TBP20  constructed the hyperparameter optimization problem through
$K$-fold cross-validation as a bilevel optimization problem with LASSO
regression and an $\ell_{1}$-norm support-vector machine (SVM) in the lower-
level problem. They used parametric programming theory to reformulate the
bilevel optimization problem as a single level problem which is called the
bilevel and parametric optimization approach to hyperparameter optimization
(HY-POP). Similar to that in TBP20 , the authors only provided the numerical
experiments to show the efficiency of HY-POP without any theoretical analysis.
KP13  considered bilevel optimization problems for variational image denoising
models, where the upper-level problem is smooth while the lower-level problem
is the $\ell_{p}$ regularizer with $p=\frac{1}{2},1,2$. They proposed
semismooth Newton method for solving the bilevel optimization problem
including the $\ell_{2}$-norm and the $\ell_{1}$-norm. Especially, they only
provided numerical experiments for the $\ell_{\frac{1}{2}}$-norm and leave the
theoretical analysis for nonconvex $\ell_{\frac{1}{2}}$-norm to future work.
Nevertheless, they showed that the $\ell_{\frac{1}{2}}$-norm has better
denoising performance than the $\ell_{1}$-norm. Recently, OTKW21  considered
the bilevel program (3) with the function
$R_{1}(\omega):=\sum_{i=1}^{n}\lvert\omega_{i}\rvert^{p}$ $(0<p\leq 1)$ (i.e.
the $\ell_{p}$-regularizer) by employing smoothing method via the twice
continuously differentiable function

$\varphi_{\mu}(\omega)=\sum_{i=1}^{n}(\omega_{i}^{2}+\mu^{2})^{\frac{p}{2}}$
(4)

as a smooth approximation of $R_{1}$. Using such a smoothing function, problem
(3) can be approximated by a smooth bilevel program, which then allows for use
of several optimization techniques which normally require differentiability.
Thus, they established the convergence analysis for nonsmooth, nonconvex
$\ell_{p}$-norm with $(0<p\leq 1)$.

In this paper, we consider a general framework for constructing smoothing
functions for $R_{1}$, where the associated $\psi$ is any function that
satisfies Assumption (A). On the one hand, this work is an important extension
of OTKW21  where the specific smoothing function used for $R_{1}$ is given by
(4) and the function $\psi$ considered is $\psi(t)=t$. However, while some
results in OTKW21  can possibly be extended to the general case, i.e. when
other smoothing functions $\varphi_{\mu}$ for $R_{1}$ are used, several
arguments and fundamental results in OTKW21  are specific to the chosen
approximation (4). For instance, (OTKW21, , Lemma 7) provides an inequality
which describes the relationship between the rates of decay to zero of the
smoothing parameter $\mu$ and the sequence of iterates $\\{\omega^{k}\\}$
produced by the smoothing algorithm. However, this can only be verified for
the smoothing function (4), and it is worthwhile to note that this proposition
is a key result to establish the convergence results. When a general smoothing
function is used for $R_{1}$, the analysis is more difficult and requires some
subtle arguments. Therefore, important contributions of this paper involve the
considerably wider class of regularizers for the hyperparameter optimization
problem, as well as the unified convergence analysis of the smoothing
algorithm for the said model under less restrictive constraint qualifications.

The smoothing framework we consider is based on CM96 , where a density
function is used to generate a smooth approximation of the plus function
$x_{+}:=\max\\{0,x\\}$, and in turn a smoothing function for $R_{1}$ (see
Section 2.3). Indeed, the theoretical novelty of our work lies on our
strategies to establish the convergence results using only our knowledge of
the density function used to obtain the smooth approximation of $R_{1}$. That
is, we do not rely on the specific formula of the smoothing function as in
OTKW21 , thus making our analysis generalized. Despite the difficulty in
accomplishing this task, we successfully identified important properties of
density functions that will suffice to prove our main results, particularly
the convergence of the smoothing algorithm to a candidate solution of (3). Of
course, the main results of the proposed smoothing framework in this present
work subsume those presented in OTKW21 . However, we emphasize that some
results in OTKW21  cannot be simply extended to the general case since in the
said work, the formula of the specific smoothing function chosen was fully
exploited.

From a practical point of view, the choice of smoothing functions is critical
in achieving successful simulations with fast convergence rates. We compare
the numerical performance of six smoothing functions generated via Chen and
Mangasarian’s method CM96  to determine which function is more suitable for
the smoothing approach.

This paper is organized as follows: In Section 2, we review some fundamental
concepts in analysis and a brief introduction about first-order necessary
optimality condition for lower-level problem. We also review the method
proposed in CM96  to construct smoothing functions of the plus function by
means of density functions. This will serve as our basis to construct
smoothing functions for $R_{1}(\omega)$, and our theoretical analysis will all
be dependent on the density function. In Section 3, the approximation of
lower-level problem and the convergence analysis are proposed for the bilevel
nonsmooth optimization problem (3). In Section 4, we list six smoothing
functions of $R_{1}(\omega)$ and compare the applicability and efficiency of
these smoothing functions in solving (3).

Throughout this paper, we denote the vector $\omega\in\Re^{n}$ by
$\omega=(\omega_{1},\dots,\omega_{n})^{T}$, and we let
$\lvert\omega\rvert:=(\lvert\omega_{1}\rvert,\dots,\lvert\omega_{n}\rvert)^{T}$,
and
$\lvert\omega\rvert^{p}:=(\lvert\omega_{1}\rvert^{p},\dots,\lvert\omega_{p}\rvert^{p})^{T}$.
We define the ${\rm sgn}$ function as sgn$(t)=1$ $(t>0)$, sgn$(t)$ = 0
$(t=0)$, and sgn$(t)=-1$ $(t<0)$. For a differentiable function
$f:\Re^{n}\to\Re$, we denote the gradient function of $f$ by $\nabla f$ with
$\nabla f(\omega):=(\frac{\partial
f(\omega)}{\partial\omega_{1}},\dots,\frac{\partial
f(\omega)}{\partial\omega_{n}})^{T}\in\Re^{n}$ and if $f$ is twice
differentiable, we denote the Hessian of $f$ by $\nabla^{2}f$ with
$\nabla^{2}f(\omega):=\left(\frac{\partial^{2}f(\omega)}{\partial\omega_{i}\partial\omega_{j}}\right)_{1\leq
i,j\leq n}\in\Re^{n\times n}$.

##  2 Preliminaries

We review some important concepts in nonsmooth analysis. We also present the
corresponding scaled bilevel Karush-Kuhn-Tucker conditions of our main problem
(3). Finally, we review the method of Chen and Mangasarian to construct
smoothing functions for the plus function and discuss how to use this to
obtain a smoothing function for the absolute value function.

###  2.1 Some concepts in analysis

The following facts can be found in the books of Rockafellar and Wet RW09 .

######  Definition 2.1.

(RW09, , Definition 8.3) Let $f:\Re^{n}\to\Re\cup\\{\infty\\}$ be a proper
function. For vectors $v\in\Re^{n}$ and $\bar{x}\in\Re^{n}$, one say that

  * 1.

$v$ is a regular subgradient of $f$ at $\bar{x}$, written
$v\in\hat{\partial}f(\bar{x})$, if

$f(x)\geq f(\bar{x})+v^{T}(x-\bar{x})+o(\|x-\bar{x}\|).$

  * 2.

$v$ is a general subgradient of $f$ at $\bar{x}$, written $v\in\partial
f(\bar{x})$, if there are sequences $\\{x^{\nu}\\}\subseteq\Re^{n}$ and
$\\{v^{\nu}\\}\subseteq\Re^{n}$ such that

$\lim\limits_{\nu\to\infty}x^{\nu}=\bar{x}\;\mbox{and}\;v^{\nu}\in\hat{\partial}f(x^{\nu})\;\mbox{with}\;\lim\limits_{\nu\to\infty}v^{\nu}=v.$

Note that a regular subgradient of $f$ at $\bar{x}$ is also called a Fréchet
subgradient of $f$ at $\bar{x}$ (see in K03 ). Moreover, if $f$ is a proper
and convex function, the regular subgradient of $f$ coincides the subgradient
of $f$ in the sense of convex analysis (see in (RW09, , Proposition 8.12)).

######  Proposition 2.1.

(RW09, , Theorem 8.6) For a function $f:\Re^{n}\to\Re\cup\\{\infty\\}$ and a
point $\bar{x}$ where $f$ is finite, the subgradient sets
$\hat{\partial}f(\bar{x})$ and $\partial f(\bar{x})$ are closed, with
$\hat{\partial}f(\bar{x})$ convex and $\hat{\partial}f(\bar{x})\subset\partial
f(\bar{x})$.

######  Proposition 2.2.

(RW09, , Theorem 10.1) If a proper function $f:\Re^{n}\to\Re\cup\\{\infty\\}$
has a local minimum at $\bar{x}$, then
$0\in\hat{\partial}f(\bar{x})\subset\partial f(\bar{x})$.

###  2.2 First-order optimality condition for the lower-level problem

Using Proposition 2.2, the first-order optimality condition for the lower-
level problem (2) is given by

$0\in\partial_{\omega}(G(\omega^{*},\bar{\lambda})+\lambda_{1}R_{1}(\omega^{*})),$

where
$\partial_{\omega}(G(\omega^{*},\bar{\lambda})+\lambda_{1}R_{1}(\omega^{*}))$
is the general subgradient with respect to $\omega$ of
$G(\omega,\bar{\lambda})+\lambda_{1}R_{1}(\omega)$ at $\omega^{*}$. Then
problem (3) can be transformed into the following one-level problem

$\begin{array}[]{cc}\min\limits_{\omega,\lambda}&f(\omega)\\\ {\rm
s.t}&0\in\partial_{\omega}(G(\omega,\bar{\lambda})+\lambda_{1}R_{1}(\omega))\\\
&\lambda\geq 0.\end{array}$ (5)

To obtain candidate solutions of problem (5), we employ the scaled first-order
necessary condition for a non-Lipschitz continuous function which was proposed
by Chen, Xu and Ye CXY10 . Since the function
$G(\omega,\bar{\lambda})+\lambda_{1}R_{1}(\omega)$ $(0<p<1)$ is non-Lipschitz
continuous function, we define the scaled first-order necessary condition for
the lower-level problem (2) as follows.

######  Definition 2.2.

We say that $\omega^{*}$ satisfies the scaled first-order necessary condition
of (2) if

$W_{*}\nabla_{\omega}G(\omega^{*},\bar{\lambda})+p\lambda_{1}\lvert
W_{*}\rvert^{p}\psi^{\prime}(\lvert\omega^{*}\rvert^{p})=0,$

where $W_{*}:={\rm diag}(\omega^{*})$, $\lvert W_{*}\rvert^{p}:={\rm
diag}(\lvert\omega^{*}\rvert^{p})$, and

$\psi^{\prime}(\lvert\omega^{*}\rvert^{p}):=(\psi^{\prime}(\lvert\omega_{1}^{*}\rvert^{p}),\psi^{\prime}(\lvert\omega_{2}^{*}\rvert^{p}),\dots,\psi^{\prime}(\lvert\omega_{n}^{*}\rvert^{p}))^{T}.$

Using this scaled first-order necessary condition, we obtain the following
one-level problem.

$\begin{array}[]{cc}\min\limits_{\omega,\lambda}&f(\omega)\\\ {\rm
s.t}&W\nabla_{\omega}G(\omega,\bar{\lambda})+p\lambda_{1}\lvert
W\rvert^{p}\psi^{\prime}(\lvert\omega\rvert^{p})=0\\\ &\lambda\geq
0,\end{array}$ (6)

where $W:={\rm diag}(\omega)$, $\lvert W\rvert^{p}:={\rm
diag}(\lvert\omega\rvert^{p})$, and

$\psi^{\prime}(\lvert\omega\rvert^{p}):=(\psi^{\prime}(\lvert\omega_{1}\rvert^{p}),\psi^{\prime}(\lvert\omega_{2}\rvert^{p}),\dots,\psi^{\prime}(\lvert\omega_{n}\rvert^{p}))^{T}.$

Notice that this problem is different from (5). The following result shows the
relationship between the feasible regions of problems (5) and (6) which can
proved in the same manner as (OTKW21, , Lemma 3).

######  Proposition 2.3.

For $\omega\in\Re^{n}$ and $\lambda\in\Re_{+}^{r}$, if
$0\in\partial_{\omega}(G(\omega,\bar{\lambda})+\lambda_{1}R_{1}(\omega))$,
then $W\nabla_{\omega}G(\omega,\bar{\lambda})+p\lambda_{1}\lvert
W\rvert^{p}\psi^{\prime}(\lvert\omega\rvert^{p})=0$. In particular, when
$p<1$, the converse is also true.

From the above proposition, we observe that the feasible region of (6) is
larger than the feasible region of (5) when $p=1$. Moreover, the feasible
regions of problems (5) and (6) are identical when $p<1$.

Now, we state a new optimality condition for problem (3) based on the scaled
first-order necessary condition, namely scaled bilevel KKT (SB-KKT) conditions
as the following definition.

######  Definition 2.3.

We say that $(\omega^{*},\lambda^{*})\in\Re^{n}\times\Re^{r}$ is a scaled
bilevel Karush-Kuhn-Tucker (SB-KKT) point for problem (3) if there exists a
pair of vectors $(\zeta^{*},\eta^{*})\in\Re^{n}\times\Re^{r}$ such that

$\displaystyle W^{2}_{*}\nabla
f(\omega^{*})+H(\omega^{*},\lambda^{*})\zeta^{*}=0,$ (7) $\displaystyle
W_{*}\nabla_{\omega}G(\omega^{*},\bar{\lambda}^{*})+p\lambda^{*}_{1}\lvert
W_{*}\rvert^{p}\psi^{\prime}(\lvert\omega^{*}\rvert^{p})=0,$ (8)
$\displaystyle p\sum_{j\notin I(\omega^{*})}{\rm
sgn}(\omega^{*}_{j})\lvert\omega_{j}^{*}\rvert^{p-1}\psi^{\prime}(\lvert\omega_{j}^{*}\rvert^{p})\zeta^{*}_{j}=\eta^{*}_{1},$
(9) $\displaystyle\zeta^{*}_{j}=0\ (j\in I(\omega^{*})),$ (10)
$\displaystyle\nabla R_{j}(\omega^{*})^{T}\zeta-\eta^{*}_{j}=0\
(j=2,3,\dots,r),$ (11) $\displaystyle
0\leq\lambda^{*},0\leq\eta^{*},(\lambda^{*})^{T}\eta^{*}=0,$ (12)

where $W_{*}:={\rm diag}(\omega^{*})$,
$I(\omega):=\\{j\in\\{1,2,\dots,n\\}~{}\lvert~{}\omega_{j}=0\\}$. Here, we
write

$H(\omega,\lambda)=W^{2}\nabla^{2}_{\omega\omega}G(\omega,\bar{\lambda})+\lambda_{1}p(p-1){\rm
diag}(\lvert
W\rvert^{p}\psi^{\prime}(\lvert\omega\rvert^{p}))+\lambda_{1}p^{2}{\rm
diag}(\lvert W\rvert^{2p}\psi^{\prime\prime}(\lvert\omega\rvert^{p}))$

with $W:={\rm diag}(\omega)$, $\lvert W\rvert^{p}:={\rm
diag}(\lvert\omega\rvert^{p})$, and $\lvert W\rvert^{2p}:={\rm
diag}(\lvert\omega\rvert^{2p})$ for $\omega\in\Re^{n}$ and
$\lambda\in\Re^{r}$.

The SB-KKT conditions are necessary optimality conditions for problem (6) as
asserted by the following results. The proof of this proposition is
essentially similar to (OTKW21, , Theorems 2,4). Hence, we only state the
result and omit the detailed proof.

######  Proposition 2.4.

Let $(\omega^{*},\lambda^{*})\in\Re^{n}\times\Re^{r}$ be a local optimum of
(6). Then, $(\omega^{*},\lambda^{*})$ together with some pair of vectors
$(\zeta^{*},\eta^{*})\in\Re^{n}\times\Re^{r}$ satisfies the SB-KKT conditions
(7)-(12) under an appropriate constraint qualification concerning the
constraints $\frac{\partial G(\omega,\bar{\lambda})}{\partial\omega_{j}}+p{\rm
sgn}(\omega_{j})\lambda_{1}\lvert\omega_{j}\rvert^{p-1}\psi^{\prime}(\lvert\omega_{j}\rvert^{p})=0\;(j\notin
I(\omega^{*}))$, $\omega_{j}=0\;(j\in I(\omega^{*}))$, and $\lambda\geq 0$.

######  Proposition 2.5.

Let $p<1$ and $(\omega^{*},\lambda^{*})\in\Re^{n}\times\Re^{r}$ be a local
optimum of (5). Then, $(\omega^{*},\lambda^{*})$ together with some pair of
vectors $(\zeta^{*},\eta^{*})\in\Re^{n}\times\Re^{r}$ satisfies the SB-KKT
conditions (7)-(12) under an appropriate constraint qualification concerning
the constraints $\frac{\partial
G(\omega,\bar{\lambda})}{\partial\omega_{j}}+p{\rm
sgn}(\omega_{j})\lambda_{1}\lvert\omega_{j}\rvert^{p-1}\psi^{\prime}(\lvert\omega_{j}\rvert^{p})=0\;(j\notin
I(\omega^{*}))$, $\omega_{j}=0\;(j\in I(\omega^{*}))$, and $\lambda\geq 0$.

###  2.3 Smoothing functions of $\lvert x\rvert$ via density functions

We next recall the general definition of a smoothing function.

######  Definition 2.4.

(C12, , Definition 1) Let $h:\Re^{n}\to\Re$ be a continuous function. We say
that $\phi:\Re_{++}\times\Re^{n}\to\Re$ is a smoothing function of $h$ if it
satisfies the following:

(i)

$\phi(\mu,\cdot)$ is continuously differentiable for any $\mu>0$;

(ii)

$\displaystyle\lim_{w\to z,\mu\downarrow 0}\phi(\mu,w)=h(z)$ for any
$z\in\Re^{n}$.

We now review the construction way of smoothing functions of the plus function
$(x)_{+}=\max\\{x,0\\}$ for $x\in\Re$ which was studied by Chen and
Mangasarian CM96 . It is well-known that
$(x)_{+}=\int^{x}_{-\infty}\sigma(t)dt$, where $\sigma$ is the step function
given by

$\sigma(x)=\left\\{\begin{array}[]{ccc}1&{\rm if}&x>0,\\\ 0&{\rm if}&x\leq
0.\end{array}\right.$

On the other hand, the step function $\sigma$ can be described as
$\sigma(x)=\int^{x}_{-\infty}\delta(t)dt$, where $\delta$ is the Dirac delta
function which satisfies the following properties

$\delta(x)=\left\\{\begin{array}[]{ccc}0&{\rm if}&x\neq 0\\\ \infty&{\rm
if}&x=0\\\ \end{array}\right.\ \mbox{and}\
\int^{\infty}_{-\infty}\delta(x)dx=1.$

Clearly, the plus function can be obtained by integrating twice the Dirac
delta function. To obtain a smoothing function of the plus function, we start
by considering a probability density function, which will serve as a smoothing
function of the Dirac delta function. More precisely, let $\rho:\Re\to\Re_{+}$
be a piecewise continuous function, called a density or kernel function, with
a finite number of pieces satisfying the following properties

$\rho(x)\geq 0\ \mbox{and}\ \int_{-\infty}^{+\infty}\rho(x)dx=1.$

Define

$\hat{t}(\mu,x):=\frac{1}{\mu}\rho\left(\frac{x}{\mu}\right),$

where $\mu$ is a positive parameter. When $\mu$ goes to 0, the limit of
$\hat{t}(\mu,x)$ is the Dirac delta function $\delta(x)$. A smooth
approximation of $\sigma$ can then be formed as

$\hat{s}(\mu,x):=\int_{-\infty}^{x}\hat{t}(\mu,t)dt\approx\sigma(x),$

and finally, this yields a smoothing function $\hat{\phi}$ for the plus
function given by

$\hat{\phi}(\mu,x):=\int_{-\infty}^{x}\hat{s}(\mu,t)dt\approx(x)_{+}.$

The above-mentioned relationship among the functions is summarized as follows:

$\displaystyle\begin{array}[]{ccccccc}&&\delta(x)&\xrightarrow{\int_{-\infty}^{x}\delta(t)dt}&\
\sigma(x)&\xrightarrow{\int_{-\infty}^{x}\sigma(t)dt}&\ (x)_{+}\vspace{1em}\\\
&&\rotatebox{90.0}{$\thickapprox$}&&\rotatebox{90.0}{$\thickapprox$}&&\rotatebox{90.0}{$\thickapprox$}\\\
\rho(x)&\xrightarrow{\frac{1}{\mu}\rho\left(\frac{x}{\mu}\right)}&\
\hat{t}(\mu,x)&\xrightarrow{\int_{-\infty}^{x}\hat{t}(\mu,t)dt}&\
\hat{s}(\mu,x)&\xrightarrow{\int_{-\infty}^{x}\hat{s}(\mu,t)dt}&\
\hat{\phi}(\mu,x)\end{array}$

In addition, assume that the function $\rho:\Re\to\Re_{+}$ satisfies

$\rho(x)=\rho(-x)\ \mbox{and}\ \kappa:=\int_{-\infty}^{+\infty}\lvert
x\rvert\rho(x)dx<+\infty,$ (13)

implying that, for any $\mu>0$,

$\hat{t}(\mu,x)=\hat{t}(\mu,-x)\ \mbox{and}\ \int_{-\infty}^{+\infty}\lvert
x\rvert\hat{t}(\mu,x)dx<+\infty.$

Therefore, the smooth approximation $\hat{\phi}$ of $(x)_{+}$ can be described
by convolution C12 ; CM96 , i.e.,

$\hat{\phi}(\mu,x)=\int_{-\infty}^{+\infty}(x-t)_{+}\hat{t}(\mu,t)dt=\int_{-\infty}^{x}(x-t)\hat{t}(\mu,t)dt.$
(14)

The following proposition shows the properties of $\hat{\phi}(\mu,x)$, whose
proof can be found in (CM96, , Proposition 2.2).

######  Proposition 2.6.

(CM96, , Proposition 2.2) Suppose that $\hat{\phi}(\mu,x)$ is defined by (14).
Then, for a fixed $\mu>0$, $\hat{\phi}(\mu,x)$ possesses the following
properties.

(a)

$\hat{\phi}(\mu,x)$ is continuously differentiable.

(b)

$-D_{2}\mu\leq\hat{\phi}(\mu,x)-(x)_{+}\leq D_{1}\mu$, where

$D_{1}:=\int_{-\infty}^{0}\lvert x\rvert\rho(x)dx,\enskip\mbox{and}\enskip
D_{2}:=\max\left\\{\int_{-\infty}^{+\infty}x\rho(x)dx,0\right\\}.$

(c)

$\hat{\phi}^{\prime}(\mu,x)$ is bounded satisfying
$0\leq\hat{\phi}^{\prime}(\mu,x)\leq 1$, where $\hat{\phi}^{\prime}(\mu,x)$ is
the first derivative of $\hat{\phi}(\mu,x)$ with respect to variable $x$.

By Definition 2.4 and Proposition 2.6, we see that $\hat{\phi}(\mu,x)$ is a
smoothing function of the plus function $(x)_{+}$. In fact, from (13) and
Proposition 2.6(b), we have that

$0\leq\hat{\phi}(\mu,x)-(x)_{+}\leq\frac{1}{2}\kappa\mu$

which implies that for any $x\in\Re$,

$\lim_{x^{k}\to x,\mu_{k}\downarrow 0}\hat{\phi}(\mu_{k},x^{k})=(x)_{+}.$ (15)

Meanwhile, since $\hat{\phi}^{\prime}(\mu,x)=\int_{-\infty}^{x/\mu}\rho(t)dt$,
then

$\lim_{x^{k}\to x,\mu_{k}\downarrow
0}\hat{\phi}^{\prime}(\mu_{k},x^{k})=\begin{cases}1&\text{if}~{}x>0\\\
0&\text{if}~{}x<0.\end{cases}$ (16)

Although the above limit may not exist when $x=0$, we have from Proposition
2.6(c) that subsequential limits exist and belong to the interval $[0,1]$.

To obtain a smoothing function for the absolute value function, observe that
$\lvert x\rvert=(x)_{+}+(-x)_{+}$. Then we achieve a smoothing function of
$\lvert x\rvert$ as follows:

$\phi(\mu,x):=\hat{\phi}(\mu,x)+\hat{\phi}(\mu,-x)=\int_{-\infty}^{+\infty}\lvert
x-t\rvert\hat{t}(\mu,t)dt.$ (17)

In a manner similar to Proposition 2.6, we also have the following properties
for $\phi(\mu,x)$.

######  Proposition 2.7.

Suppose that $\phi(\mu,x)$ is defined as in (17). Then, for a fixed $\mu>0$,
we have

(a)

$\phi(\mu,x)$ is continuously differentiable.

(b)

$0\leq\phi(\mu,x)-\lvert x\rvert\leq\kappa\mu$ with the constant $\kappa>0$
defined in (13).

(c)

$\phi^{\prime}(\mu,x)$ is bounded satisfying $-1\leq\phi^{\prime}(\mu,x)\leq
1$.

From Proposition 2.7, we obtain results analogous to (15) and (16). That is,

$\lim_{x^{k}\to x,\mu_{k}\downarrow 0}\phi(\mu_{k},x^{k})=\lvert
x\rvert,\quad\forall x\in\Re$

and

$\lim_{x^{k}\to x,\mu_{k}\downarrow 0}\phi^{\prime}(\mu_{k},x^{k})={\rm
sgn}(x)\quad\forall x\neq 0.$ (18)

Likewise, we have from Proposition 2.7(c) that subsequential limits of the
sequence on the left-hand side of (18) exist and belong to $[-1,1]$.

There are various choices of the density function $\rho$ and corresponding
smoothing function $\phi$. Many instances will be presented in the section of
numerical experiments.

##  3 Proposed algorithm and convergence results

In this section, we describe our smoothing algorithm for (3) and present our
convergence results.

###  3.1 Smoothing approach and the algorithm

One main source of difficulty in solving the bilevel program (3) is the
nonsmooth, nonconvex and possibly non-Lipschitz component
$R_{1}(\omega)=\sum_{i=1}^{n}\psi(\lvert\omega_{i}\rvert^{p})$, where
$p\in(0,1]$. To overcome this difficulty, we apply the smoothing technique to
$R_{1}$ with the smoothing function $\phi$ defined in the previous section,
yielding the following smoothed function:

$\varphi_{\mu}(\omega):=\sum_{j=1}^{n}\psi\left([\phi(\mu,\omega_{j})]^{p}\right).$
(19)

Then, as in OTKW21 , we consider problem (3) with $\varphi_{\mu}$ in place of
$R_{1}$, and further replace the obtained smoothed lower-level problem with
its first-order condition. Hence, the following problem is obtained:

$\begin{array}[]{cc}\min\limits_{\omega,\lambda}&f(\omega)\\\ {\rm
s.t}&\nabla_{\omega}G(\omega,\bar{\lambda})+\lambda_{1}\nabla\varphi_{\mu}(\omega)=0\\\
&\lambda\geq 0.\end{array}$ (20)

Next, let us suppose that $\varphi_{\mu}$ is twice continuously differentiable
from this moment and recall Assumption (A) together with our differentiability
assumptions on $R_{i}$ ($i=2,\dots,r$) and $g$. These properties enable us to
consider the KKT conditions.111Without the $\mathcal{C}^{2}$ property of
$\varphi_{\mu}$, the KKT conditions cannot be taken because the constraint
function
$\nabla_{\omega}G(\omega,\bar{\lambda})+\lambda_{1}\nabla\varphi_{\mu}(\omega)$
is not necessarily smooth. By virtue of this fact, we can find candidate
solutions to (20) by looking at its KKT points.

In fact, it is sufficient to obtain approximate KKT points as introduced in
OTKW21 : Given a parameter $\hat{\varepsilon}>0$, we define an
$\hat{\varepsilon}$-approximate KKT point for problem (20) as follows: We say
that
$\\{(\omega,\lambda,\zeta,\eta)\\}\subseteq\Re^{n}\times\Re^{r}\times\Re^{n}\times\Re^{r}$
is an $\hat{\varepsilon}$-approximate KKT point for (20) if there exists a
vector
$(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4},\varepsilon_{5})\in\Re^{n}\times\Re\times\Re^{r-1}\times\Re^{n}\times\Re$
such that

$\displaystyle\nabla
f(\omega)+(\nabla^{2}_{\omega\omega}G(\omega,\bar{\lambda})+\lambda_{1}\nabla^{2}\varphi_{\mu}(\omega))\zeta=\varepsilon_{1},$
(21)
$\displaystyle\nabla\varphi_{\mu}(\omega)^{T}\zeta-\eta_{1}=\varepsilon_{2},$
(22) $\displaystyle\nabla
R_{j}(\omega)^{T}\zeta-\eta_{j}=(\varepsilon_{3})_{j}\ (j=2,3,\dots,r),$ (23)
$\displaystyle\nabla_{\omega}G(\omega,\bar{\lambda})+\lambda_{1}\nabla\varphi_{\mu}(\omega)=\varepsilon_{4},$
(24) $\displaystyle 0\leq\lambda,\quad
0\leq\eta,\quad\lambda^{T}\eta=\varepsilon_{5},$ (25)

and

$\|(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4},\varepsilon_{5})\|\leq\hat{\varepsilon},$

where $\nabla^{2}_{\omega\omega}G(\omega,\bar{\lambda})$ is the Hessian of $G$
with respect to $\omega$. Note that when $\hat{\varepsilon}=0$, an
$\hat{\varepsilon}$-approximate KKT point is identical to a KKT point.

Now, by iteratively computing an $\hat{\varepsilon}$-approximate KKT point
while decreasing $\hat{\varepsilon}$ and the smoothing parameter $\mu$, we
obtain the following smoothing algorithm:

Algorithm 1 (A Smoothing Method for Nonsmooth Bilevel Optimization)

Step 0

Choose $\mu_{0}\neq 0$, $\beta_{1},\beta_{2}\in(0,1)$ and
$\hat{\varepsilon}_{0}\geq 0$. Set $k:=0$.

Step 1

Find an $\hat{\varepsilon}_{k}$-approximate KKT point
$\\{(\omega^{k+1},\lambda^{k+1},\zeta^{k+1},\eta^{k+1})\\}$ for problem (20)
with $\mu=\mu_{k}$.

Step 2

Set $\mu_{k+1}=\beta_{1}\mu_{k}$,
$\hat{\varepsilon}_{k+1}=\beta_{2}\hat{\varepsilon}_{k}$ and $k:=k+1$.

Algorithm 1 is quite similar to the one proposed by Okuno et al OTKW21 .
However, whereas Okuno et al supposed to employ only
$\sum_{i=1}^{n}(\omega_{i}^{2}+\mu^{2})^{\frac{p}{2}}$ as the smoothing
function $\varphi_{\mu}$, it enjoys much more freedom in choices.

In order to establish the global convergence of Algorithm 1, we will require
some more properties of the density function $\rho$ used for constructing the
smoothing function $\phi$, which also accomplish the
$\mathcal{C}^{2}$-property of $\varphi_{\mu}$ supposed above.

###  3.2 Assumptions for convergence analysis

####  3.2.1 Assumptions on iterates

We assume that at every iteration, an $\hat{\varepsilon}_{k}$-approximate KKT
point can always be computed henceforth. We also make the following technical
assumptions about the sequence generated by Algorithm 1:

Assumption B. Let
$\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}\subseteq\Re^{n}\times\Re^{r}\times\Re^{n}\times\Re^{r}$
be the sequence produced by Algorithm 1. Then, the following properties hold,
where $(\omega^{*},\lambda^{*})$ denotes an arbitrary accumulation point of
the sequence $\\{(\omega^{k},\lambda^{k})\\}$:

  * (B1)

$\liminf_{k\to\infty}\lambda_{1}^{k}>0$.

  * (B2)

The sequence $\\{(\omega^{k},\lambda^{k})\\}$ is bounded.

  * (B3)

Let $p=1$ and $\lambda^{*}=(\lambda_{1}^{*},\bar{\lambda}^{*})$. Assume that
$\lambda^{*}_{1}\psi^{\prime}(0)\neq\lvert\frac{\partial
G(\omega^{*},\bar{\lambda}^{*})}{\partial\omega_{j}}\rvert$ for any $j\in
I(\omega^{*})$, where $I(\omega^{*})$ is described in Definition 2.3.

  * (B4)

Denote

$I(\lambda^{*}):=\\{i\in\\{1,2,\dots,r\\}\ \lvert\ \lambda^{*}_{i}=0\\}$

and

$\Phi_{j}(\omega,\lambda):=\frac{\partial
G(\omega,\bar{\lambda})}{\partial\omega_{j}}+p{\rm
sgn}(\omega_{j})\lambda_{1}\lvert\omega_{j}\rvert^{p-1}\psi^{\prime}(\lvert\omega_{j}\rvert^{p})\;(j\notin
I(\omega^{*})).$

Then, the Mangasarian-Fromovitz constraint qualification (MFCQ) holds at
$(\omega,\lambda)=(\omega^{*},\lambda^{*})$ for the constraints
$\Psi_{j}(\omega,\lambda)=0$ for all $j=1,\dots,n$ and $\lambda\geq 0$, where

$\Psi_{j}(\omega,\lambda):=\left\\{\begin{array}[]{lc}\Phi_{j}(\omega,\lambda)&{\rm
if}\,j\notin I(\omega^{*}),\\\ \omega_{j}&{\rm if}\,j\in I(\omega^{*}).\\\
\end{array}\right.$

That is,
$\\{\nabla_{(\omega,\lambda)}\Psi_{j}(\omega^{*},\lambda^{*})\\}_{j=1}^{n}$ is
linearly independent and there exists $\bar{d}\in\Re^{n+r}$ such that

$\displaystyle\nabla_{(\omega,\lambda)}\Psi_{j}(\omega^{*},\lambda^{*})^{T}\bar{d}=0\quad\forall
j=1,\dots,n,$ (26)
$\displaystyle(\nabla_{(\omega,\lambda)}\lambda_{i}\lvert_{(\omega,\lambda)=(\omega^{*},\lambda^{*})})^{T}\bar{d}>0\qquad\forall
i\in I(\lambda^{*}).$ (27)

Note that Assumption (B4) is weaker than the linearly independent constraint
qualification (LICQ) which was used in OTKW21 . The following lemma is needed
for subsequent analysis.

######  Lemma 3.1.

Suppose that $(\omega^{*},\lambda^{*})$ is an arbitrary accumulation point of
the sequence $\\{(\omega^{k},\lambda^{k})\\}$ such that Assumption (B4) holds.
Then,
$\\{\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})\\}_{j\notin
I(\omega^{*})}$ is linearly independent and there exists a vector
$d\in\Re^{n-\lvert I(\omega^{*})\rvert+r}$ such that

$\displaystyle\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})^{T}d=0\quad\forall
j\notin I(\omega^{*}),$ (28)
$\displaystyle(\nabla_{(\bar{\omega},\lambda)}\lambda_{i}\lvert_{(\omega,\lambda)=(\omega^{*},\lambda^{*})})^{T}d>0\quad\forall
i\in I(\lambda^{*}),$ (29)

where $\bar{\omega}:=(\omega_{j})_{j\in I(\omega^{*})}$.

Proof: See Appendix B.

####  3.2.2 Assumptions on density function

We have mentioned in the Introduction that the setting of all our analysis is
based on density functions. That is, we wish to prove all our convergence
results by solely looking at density functions used to induce the smoothing
functions. To this end, we must be able to identify necessary properties of a
given density function so that Algorithm 1 converges to a candidate solution
of the main problem (3). Indeed, the novelty of our work is precisely the
identification of these required properties and its application to the
convergence analysis.

We summarize necessary assumptions on the density function $\rho$ that we will
use in the next subsection.

Assumption C. Let $\rho:\Re\to\Re_{+}$ be a density function. Then, the
following properties hold:

  * (C1)

$\rho$ is symmetric, i.e. $\rho(x)=\rho(-x)$ for all $x\in\Re$.

  * (C2)

$\rho$ is continuous and nonincreasing on $[0,\infty)$.

  * (C3)

There exist positive constants $c,r>0$ such that

$2\int\limits_{0}^{S}\rho(x)~{}dx\geq 1-\frac{c}{S^{r}+c}~{}~{}\mbox{for
all}~{}S\geq 0.$

  * (C4)

If $p=1$, we have $\rho(x)>0$ for all $x\in\Re$.

Some remarks are in order: First, although Assumption (C1) was already
supposed in Section 2.3, we have restated it for later use. Under this
assumption, note that the smoothing function $\phi(\mu,x)$ is strictly
positive for all $\mu>0$ and $x\in\Re$. Indeed, we already have from
Proposition 2.7(b) that $\phi(\mu,x)>0$ for $x\neq 0$. On the other hand,
since $\rho$ is a symmetric density function by Assumption (C1), then $\rho$
is not identical to the zero function on the interval $[0,\infty)$.
Consequently, $\int_{0}^{+\infty}t\rho(t)dt>0$, which together with (14) and
(17) yields $\phi(\mu,0)>0$. Moreover, we can easily calculate the first and
second derivatives of the induced $\phi(\mu,x)$ as

$\phi^{\prime}(\mu,x)=2\mathrm{sgn}(x)\int_{0}^{\lvert
x\rvert}\frac{1}{\mu}\rho\left(\frac{t}{\mu}\right)\;dt=2\mathrm{sgn}(x)\int_{0}^{\frac{\lvert
x\rvert}{\mu}}\rho(t)\;dt,$ (30)

and

$\phi^{\prime\prime}(\mu,x)=\frac{2}{\mu}\rho\left(\frac{x}{\mu}\right),$ (31)

respectively. From equation (31) and strict positivity of $\phi(\mu,x)$, we
see that $\phi(\mu,\cdot)$ is twice continuously differentiable approximation
of $\lvert x\rvert$ by the continuity assumption in (C2).

Assumptions (C1) and (C2) will also have other important roles in the proofs
of our main result. On the other hand, the other two technical assumptions on
$\rho$, namely (C3) and (C4), are important in our subsequent analysis.
Without knowing definitively the formula for $\phi(\mu,x)$, the analysis
becomes extremely difficult. In particular, it is challenging to understand
the convergence behavior of the sequence $\left\\{S_{j}^{k}\right\\}$, where

$S_{j}^{k}:=\frac{\lvert\omega_{j}^{k}\rvert}{\mu_{k-1}},$ (32)

$\omega^{k}=(\omega_{1}^{k},\dots,\omega_{n}^{k})$ and $\mu_{k-1}$ are
generated from Algorithm 1 when $\omega_{j}^{k}\to 0$ as $k\to\infty$ for some
$j\in\\{1,\dots,n\\}$. Nonetheless, this problem can be alleviated thanks to
the simple Assumption (C3). On the other hand, Assumption (C4) will later be
important in proving the unboundedness of the sequence
$\left\\{\lvert\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\rvert\right\\}$
when $p=1$ and $\omega_{j}^{k}\to 0$ as $k\to\infty$. Interestingly, we shall
see shortly that (C4) is not needed for the case $p\in(0,1)$.

In the next subsection, we will see in great detail how these assumptions on
$\rho$ will play a central role in establishing the main convergence result.
In Section 4.1, we will provide some specific examples of density functions
satisfying Assumption (C).

###  3.3 Convergence analysis

We now prove our main result that accumulation points of the sequence
generated by Algorithm 1 are in fact SB-KKT points (see Definition 2.3), which
in turn are candidate solutions for (6).

######  Theorem 3.1.

Any accumulation point of the sequence
$\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}$ satisfies the SB-KKT
conditions (7)-(12) for the original problem (3) provided that Assumptions
(A), (B) and (C) hold.

An accumulation point of $\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}$
will be denoted by $(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*})$. To prove
Theorem 3.1, we need to show that
$(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*})$ satisfies equations (7)-(12). We
show this by proving a series of lemmas. In particular, we do the following:

(i)

Prove that $\\{S_{j}^{k}\\}_{k\in K}$ is bounded, where $S_{j}^{k}$ is given
by (32), $j\in I(\omega^{*})$, $K\subset\\{1,2,\dots,\\}$ such that
$(\omega^{k},\lambda^{k})\to(\omega^{*},\lambda^{*})$ as $k\in K\to\infty$ and
$\\{(\omega^{k},\lambda^{k})\\}$ is generated by Algorithm 1. Such an index
set $K$ exists by Assumption (B2);

(ii)

Compute the limit of the sequence
$\left\\{(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\right\\}_{k\in K}$,
where $j\in I(\omega^{*})$ and the index set $K$ is as described in (i);

(iii)

Compute the limit of the sequences
$\left\\{(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\right\\}$ and
$\left\\{(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\right\\}$, where
$j\notin I(\omega^{*})$;

(iv)

Prove that the sequence $\\{(\zeta^{k},\eta^{k})\\}$ generated by Algorithm 1
is bounded, and thus, accumulation points of the sequence of iterates exist;

(v)

Prove that any accumulation point
$(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*})$ of the sequence
$\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}$ generated by Algorithm 1
satisfies equations (9) and (10); and finally,

(vi)

Compute the limits of the sequences
$\left\\{W_{k}\nabla\varphi_{\mu_{k-1}}(\omega^{k})\right\\}_{k\in K}$ and
$\left\\{W_{k}^{2}\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k})\right\\}_{k\in K}$
for an arbitrary sequence $\\{(\omega^{k},\mu_{k-1})\\}$ converging to
$(\omega^{*},0)$, where $W_{k}:={\rm diag}(\omega^{k})$.

The above objectives are formally stated and proved, respectively, in Lemma
3.2 to Lemma 3.7. We will prove these results without knowledge of the exact
formula for the smoothing function $\phi(\mu,x)$ used to construct
$\varphi_{\mu}$ given by (19), that is, only using Assumption (C) on the
density function.

In OTKW21 , the authors proved that when $\psi(t)=t$ and
$\phi(\mu,x)=\sqrt{x^{2}+\mu^{2}}$ and under Assumptions (B1)-(B3), there
exists some $\gamma>0$ such that

$\mu_{k-1}^{2}\geq\gamma\lvert\omega_{j}^{k}\rvert^{\frac{2}{2-p}}\quad(j\in
I(\omega^{*}))$ (33)

for all sufficiently large $k\in K$, where $K$ is the subsequence described in
the above (i). This result is especially important in proving results
analogous to (ii)-(iv) above. However, in order to derive inequality (33),
OTKW21  takes advantage of the specific function $\phi(\mu,x)$ chosen, which
is not the case in the present manuscript. Nevertheless, we have found out
that such a strong result is not necessarily required to prove (ii)-(iv). In
particular, it suffices to establish (i), which is indeed a weaker property.
To this end, Assumption (C3) will play a very significant role without which
the analysis becomes extremely difficult.

First of all, it can be easily calculated that the components of
$\nabla\varphi_{\mu}(\omega)\in\Re^{n}$ are given by

$(\nabla\varphi_{\mu}(\omega))_{j}=p\psi^{\prime}\left([\phi(\mu,\omega_{j})]^{p}\right)\phi^{\prime}(\mu,\omega_{j})[\phi(\mu,\omega_{j})]^{p-1},$
(34)

while $\nabla^{2}\varphi_{\mu}(\omega)\in\Re^{n\times n}$ is a diagonal matrix
whose diagonal entries are given by

$\displaystyle(\nabla^{2}\varphi_{\mu}(\omega))_{jj}$
$\displaystyle=p^{2}\psi^{\prime\prime}\left([\phi(\mu,\omega_{j})]^{p}\right)[\phi^{\prime}(\mu,\omega_{j})[\phi(\mu,\omega_{j})]^{p-1}]^{2}{}$
$\displaystyle+p(p-1)\psi^{\prime}\left([\phi(\mu,\omega_{j})]^{p}\right)[\phi^{\prime}(\mu,\omega_{j})]^{2}[\phi(\mu,\omega_{j})]^{p-2}{}$
$\displaystyle+p\psi^{\prime}\left([\phi(\mu,\omega_{j})]^{p}\right)[\phi(\mu,\omega_{j})]^{p-1}\phi^{\prime\prime}(\mu,\omega_{j}),$
(35)

for $j=1,\dots,n$. We note that as discussed in Section 3.2.2, the smoothing
function $\phi(\mu,x)$ is strictly positive under Assumption (C1) and thus,
the factor $[\phi(\mu,\omega_{j})]^{p-1}$ that appears in (34) and (35) is
real-valued. In addition, it is clear from (35) that the components of the
Hessian of $\varphi_{\mu}$ are continuous by Assumption (A), equation (31),
and Assumptions (C1)-(C2).

We now prove our first lemma which establishes property (i).

######  Lemma 3.2.

Suppose that Assumptions (B1)-(B3) and (C1)-(C3) hold. Let
$(\omega^{*},\lambda^{*})$ be an arbitrary accumulation point of the sequence
$\\{(\omega^{k},\lambda^{k})\\}$ generated by Algorithm 1, and let
$\\{(\omega^{k},\lambda^{k})\\}_{k\in K}$ be an arbitrary subsequence
converging to $(\omega^{*},\lambda^{*})$. If $j\in I(\omega^{*})$, then
$\\{S_{j}^{k}\\}_{k\in K}$ is bounded. Moreover, $S_{j}^{k}\to 0$ as $k\in
K\to\infty$ if $p\in(0,1)$.

Proof: Denote

$F_{j}(\omega^{k},\bar{\lambda}^{k}):=\frac{\partial
G(\omega^{k},\bar{\lambda}^{k})}{\partial\omega_{j}}.$ (36)

From equation (24) and (34), we have

$F_{j}(\omega^{k},\bar{\lambda}^{k})+p\lambda_{1}^{k}\psi^{\prime}\left([\phi(\mu_{k-1},\omega_{j}^{k})]^{p}\right)\phi^{\prime}(\mu_{k-1},\omega_{j}^{k})[\phi(\mu_{k-1},\omega_{j}^{k})]^{p-1}=(\varepsilon_{4}^{k-1})_{j}.$
(37)

Case 1. Suppose that $p=1$. Then

$F_{j}(\omega^{k},\bar{\lambda}^{k})+\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))\phi^{\prime}(\mu_{k-1},\omega_{j}^{k})=(\varepsilon_{4}^{k-1})_{j}.$
(38)

Rearranging the terms and using Assumptions (A) and (C3), there are constants
$c,r>0$ such that for sufficiently large $k$,

$\frac{\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}{\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}=\lvert\phi^{\prime}(\mu_{k-1},\omega_{j}^{k})\rvert=2\int_{0}^{S_{j}^{k}}\rho(s)~{}ds\geq
1-\frac{c}{(S_{j}^{k})^{r}+c},$

where the second equality holds by (30), and $\lambda_{1}^{k}>0$ for
sufficiently large $k$ by Assumption (B1). Consequently, we get

$\begin{array}[]{rl}\displaystyle\frac{c}{(S_{j}^{k})^{r}+c}&\displaystyle\geq
1-\frac{\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}{\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}\\\
\displaystyle\hfil&\displaystyle=\frac{\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))-\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}{\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}\end{array}$

Note that
$\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))-\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert>0$
for all large $k$ by Assumption (B3), Proposition 2.7(c) and (38). Then

$\displaystyle 0\leq\frac{(S_{j}^{k})^{r}}{c}\leq\frac{(S_{j}^{k})^{r}+c}{c}$
$\displaystyle\leq\frac{\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}{\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))-\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}$
$\displaystyle\to\frac{\lambda_{1}^{*}\psi^{\prime}(0)}{\lambda_{1}^{*}\psi^{\prime}(0)-\lvert
F_{j}(\omega^{*},\bar{\lambda}^{*})\rvert}~{}\text{as}~{}k\in K\to\infty$

where the finiteness of the limit is guaranteed by Assumption (B3). Hence, it
easily follows that $\\{S_{j}^{k}\\}_{k\in K}$ is bounded.

Case 2. Now, suppose $p\in(0,1)$. From equation (37), we have

$\begin{array}[]{rl}\displaystyle\frac{\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}{p\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}=\lvert\phi^{\prime}(\mu_{k-1},\omega_{j}^{k})\rvert\cdot[\phi(\mu_{k-1},\omega_{j}^{k})]^{p-1}.\end{array}$

Using Assumptions (C1) and (C3) and by (30), we get

$\displaystyle[\phi(\mu_{k-1},\omega_{j}^{k})]^{1-p}\cdot\frac{\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}{p\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}=\lvert\phi^{\prime}(\mu_{k-1},\omega_{j}^{k})\rvert\geq
1-\frac{c}{(S_{j}^{k})^{r}+c}\geq 0.$

Meanwhile, note that since $1-p>0$ and $\omega_{j}^{k}\to 0$ as $k\in
K\to\infty$, we have

$\lim_{k\in K\to\infty}[\phi(\mu_{k-1},\omega_{j}^{k})]^{1-p}=0.$

Since $\lambda_{1}^{\ast}>0$ by Assumption (B1) and $\psi^{\prime}(0)>0$ by
Assumption (A), then

$\lim_{k\in
K\to\infty}\frac{\lvert(\varepsilon_{4}^{k-1})_{j}-F_{j}(\omega^{k},\bar{\lambda}^{k})\rvert}{p\lambda_{1}^{k}\psi^{\prime}(\phi(\mu_{k-1},\omega_{j}^{k}))}=\frac{\lvert
F_{j}(\omega^{*},\bar{\lambda}^{*})\rvert}{p\lambda_{1}^{*}\psi^{\prime}(0)}.$

Thus,

$\lim_{k\in K\to\infty}\frac{c}{(S_{j}^{k})^{r}+c}=1.$

It follows that $\lim_{k\in K\to\infty}S_{j}^{k}=0$, as desired. This
completes the proof.

We now focus on the sequence
$\left\\{\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k})\right\\}_{k\in K}$. Using
Lemma 3.2, we prove the following important result, which will later be used
to establish the boundedness of the sequence $\\{\zeta^{k},\eta^{k}\\}$
generated by Algorithm 1.

######  Lemma 3.3.

Suppose that Assumptions (B1)-(B3) and (C1)-(C4) hold. Let
$(\omega^{*},\lambda^{*})$ be an arbitrary accumulation point of the sequence
$\\{(\omega^{k},\lambda^{k})\\}$ generated by Algorithm 1 and let
$\\{(\omega^{k},\lambda^{k})\\}_{k\in K}$ be an arbitrary subsequence
converging to $(\omega^{*},\lambda^{*})$. Then,

$\lim_{k\in
K\to\infty}\lvert(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\rvert=\infty\;\mbox{for}\;j\in
I(\omega^{*}).$

Proof: We first consider the case when $p=1$. In this instance, we have from
equation (35) and (31) that

$(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}=\phi^{\prime\prime}({\mu_{k-1}},\omega_{j}^{k})=\frac{2}{\mu_{k-1}}\rho\left(\frac{\omega^{k}_{j}}{\mu_{k-1}}\right).$

By Lemma 3.2, there exists $M>0$ such that
$\frac{\lvert\omega_{j}^{k}\rvert}{\mu_{k-1}}\leq M$ for all $k\in K$. Since
$\rho$ is nonincreasing on $[0,\infty)$ by Assumption (C2), we have

$\lim_{k\in
K\to\infty}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}=\lim_{k\in
K\to\infty}\frac{2}{\mu_{k-1}}\rho\left(\frac{\omega^{k}_{j}}{\mu_{k-1}}\right)\geq\lim_{k\in
K\to\infty}\frac{2}{\mu_{k-1}}\rho(M)=\infty,$

where the rightmost equality holds since $\rho(t)>0$ on $\Re$ by Assumption
(C4). This proves the claim for $p=1$.

We now consider the case when $0<p<1$. We consider two disjoint subsets of $K$

$U^{j}_{1}:=\\{k\in K~{}\lvert~{}\omega_{j}^{k}=0\\},\quad\text{and}\quad
U^{j}_{2}:=\\{k\in K~{}\lvert~{}\omega_{j}^{k}\neq 0\\},$

and the corresponding subsequences. For $k\in U^{j}_{1}$, we get from (30)
that $\phi^{\prime}(\mu_{k-1},0)=0$. From (35), we have

$(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}=p\psi^{\prime}\left([\phi(\mu_{k-1},\omega^{k}_{j})]^{p}\right)[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j}).$
(39)

Meanwhile,

$\lim_{k\in
U^{j}_{1}\to\infty}\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})=\lim_{k\in
U^{j}_{1}\to\infty}\frac{2}{\mu_{k-1}}\rho\left(\frac{\omega^{k}_{j}}{\mu_{k-1}}\right)=\lim_{k\in
U^{j}_{1}\to\infty}\frac{2}{\mu_{k-1}}\rho(0)=\infty,$ (40)

where we note that $\rho(0)>0$ by Assumption (C1) and the definition of
density function. Moreover, it is clear that

$\lim_{k\in U^{j}_{1}\to\infty}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}=\infty.$
(41)

It follows from (39), (40), (41) and Assumption (A) that

$\lim_{k\in
U^{j}_{1}\to\infty}\lvert(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\rvert=\infty.$

For $k\in U^{j}_{2}$, we obtain

$\lim_{k\in
U^{j}_{2}\to\infty}\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})=\lim_{k\to\infty}\frac{2}{\mu_{k-1}}\rho\left(\frac{\omega^{k}_{j}}{\mu_{k-1}}\right)=\infty.$
(42)

by using equation (31) and the facts that $\lim_{k\in
U^{j}_{2}\to\infty}\frac{\lvert\omega^{k}_{j}\rvert}{\mu_{k-1}}=0$ from Lemma
3.2 and $\rho$ is continuous by Assumption (C2).

Meanwhile, $\phi^{\prime\prime}(\mu_{k-1},\omega_{j}^{k})\neq 0$ for
sufficiently large $k$ using Lemma 3.2, equation (31) and Assumption (C1).
Thus, invoking Assumption (A), we have for large $k\in U_{2}^{j}$ that

$\displaystyle\psi^{\prime\prime}\left([\phi(\mu_{k-1},\omega^{k}_{j})]^{p}\right)\frac{[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}}{\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})}{}$
$\displaystyle=\psi^{\prime\prime}\left([\phi(\mu_{k-1},\omega^{k}_{j})]^{p}\right)\frac{[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p}}{\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})\phi(\mu_{k-1},\omega^{k}_{j})}{}$
$\displaystyle\geq-\beta\frac{[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p}}{\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})\phi(\mu_{k-1},\omega^{k}_{j})}.$
(43)

Using (31) and Proposition 2.7(b), we have

$\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})\phi(\mu_{k-1},\omega^{k}_{j})\geq\frac{2}{\mu_{k-1}}\rho\left(\frac{\omega^{k}_{j}}{\mu_{k-1}}\right)\lvert\omega^{k}_{j}\rvert=2S_{j}^{k}\rho(S_{j}^{k}).$

Using this fact together with (43) and equation (30), we obtain

$\displaystyle\psi^{\prime\prime}\left([\phi(\mu_{k-1},\omega^{k}_{j})]^{p}\right)\frac{[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}}{\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})}{}$
$\displaystyle\geq-\beta\frac{2\left[\int_{0}^{S_{j}^{k}}\rho(t)dt\right]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p}}{S_{j}^{k}\rho(S_{j}^{k})}{}$
$\displaystyle\geq-\beta\frac{2S_{j}^{k}(\rho(0))^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p}}{\rho(S_{j}^{k})}~{}\to~{}0~{}\mbox{as}~{}k\in
U^{j}_{2}\to\infty,$ (44)

where the last inequality holds since
$\left[\int_{0}^{S_{j}^{k}}\rho(t)dt\right]^{2}\leq(S_{j}^{k})^{2}(\rho(0))^{2}$.
Similarly, we also have

$\displaystyle(p-1)\frac{[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}}{\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})\phi(\mu_{k-1},\omega^{k}_{j})}$
$\displaystyle\geq(p-1)\frac{2\left[\int_{0}^{S_{j}^{k}}\rho(t)dt\right]^{2}}{S_{j}^{k}\rho(S_{j}^{k})}{}$
$\displaystyle\geq(p-1)\frac{2S_{j}^{k}(\rho(0))^{2}}{\rho(S_{j}^{k})}~{}\to~{}0~{}\mbox{as}~{}k\in
U^{j}_{2}\to\infty.$ (45)

For brevity, denote $\phi_{k}:=\phi(\mu_{k-1},\omega^{k}_{j})$,
$\phi^{\prime}_{k}:=\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})$,
$\phi^{\prime\prime}_{k}:=\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})$. Then
(35) can be written as

$(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}=p\phi^{\prime\prime}_{k}\phi_{k}^{p-1}\left[p\psi^{\prime\prime}\left(\phi_{k}^{p}\right)\frac{[\phi^{\prime}_{k}]^{2}\phi_{k}^{p-1}}{\phi^{\prime\prime}_{k}}+(p-1)\psi^{\prime}\left(\phi_{k}^{p}\right)\frac{[\phi^{\prime}_{k}]^{2}}{\phi^{\prime\prime}_{k}\phi_{k}}+\psi^{\prime}\left(\phi_{k}^{p}\right)\right].$
(46)

Since $0<\psi^{\prime}(t)\leq\alpha$ by Assumption (A), we have

$\lim_{k\in
U^{j}_{2}\to\infty}\left[p\psi^{\prime\prime}\left(\phi_{k}^{p}\right)\frac{[\phi^{\prime}_{k}]^{2}\phi_{k}^{p-1}}{\phi^{\prime\prime}_{k}}+(p-1)\psi^{\prime}\left(\phi_{k}^{p}\right)\frac{[\phi^{\prime}_{k}]^{2}}{\phi^{\prime\prime}_{k}\phi_{k}}+\psi^{\prime}\left(\phi_{k}^{p}\right)\right]>0$

by using the obtained limits (44) and (45). On the other hand, it is clear
from (42) that $\phi^{\prime\prime}_{k}\phi_{k}^{p-1}\to\infty$ as $k\in
U^{j}_{2}\to\infty$. Hence, taking the limit in (46),

$\lim_{k\in
U^{j}_{2}\to\infty}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}=\infty.$

This completes the proof.

We next prove the following lemma which will be used to prove the boundedness
of the iterates $\\{(\zeta^{k},\eta^{k})\\}$. In this result, we do not
actually need $\\{\omega^{k},\mu_{k-1}\\}$ to be generated from the algorithm.

######  Lemma 3.4.

Suppose that Assumption (C1)-(C2) hold. Let
$\\{(\omega^{k},\mu_{k-1})\\}\subseteq\Re^{n}\times\Re_{++}$ be an arbitrary
sequence converging to $(\omega^{*},0)$. Then for $j\notin I(\omega^{*})$,

$\displaystyle\lim_{k\to\infty}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}$
$\displaystyle=p\mathrm{sgn}(\omega^{*}_{j})\lvert\omega^{*}_{j}\rvert^{p-1}\psi^{\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right),$
(47)
$\displaystyle\lim_{k\to\infty}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}$
$\displaystyle=p^{2}\psi^{\prime\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)\lvert\omega^{*}_{j}\rvert^{2p-2}+p(p-1)\psi^{\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)\lvert\omega^{*}_{j}\rvert^{p-2}.$
(48)

Proof: Let $j\notin I(\omega^{*})$. We have from (18) that

$\lim_{k\to\infty}\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})=\mathrm{sgn}(\omega^{*}_{j}).$

Since $\phi(\mu,x)$ is a smoothing function of $\lvert x\rvert$ that is
strictly positive by Assumption (C1), we have
$[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}\to\lvert\omega^{*}_{j}\rvert^{p-1}$ as
$k\to\infty$. Moreover, since $\psi$ is $\mathcal{C}^{2}$ by Assumption (A),
then we easily obtain (47) by letting $k\to\infty$ in equation (34).

For equation (48), we first show that
$\lim_{k\to\infty}\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})=0$, which by
(31) is equivalent to showing that

$\lim_{k\to\infty}\frac{1}{\mu_{k-1}}\rho\left(\frac{\omega^{k}_{j}}{\mu_{k-1}}\right)=\lim_{k\to\infty}\frac{1}{\mu_{k-1}}\rho\left(\frac{\lvert\omega^{k}_{j}\rvert}{\mu_{k-1}}\right)=0.$
(49)

Since $j\notin I(\omega^{*})$, then
$\lvert\omega^{*}_{j}\rvert/2<\lvert\omega^{k}_{j}\rvert$ for all $k$
sufficiently large. Thus,

$0\leq\frac{1}{\mu_{k-1}}\rho\left(\frac{\omega_{j}^{k}}{\mu_{k-1}}\right)\leq\frac{1}{\mu_{k-1}}\rho\left(\frac{\omega_{j}^{*}}{2\mu_{k-1}}\right)\to
2\delta(\omega_{j}^{*})=0\qquad\text{as}~{}k\to\infty,$

by invoking Assumptions (C1)-(C2) and the definition of the Dirac delta
function. This proves (49). Finally, taking the limit in (35) when
$k\to\infty$, we obtain

$\displaystyle\lim_{k\to\infty}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}$
$\displaystyle=p^{2}\psi^{\prime\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)[\mbox{sgn}(\omega^{*}_{j})\lvert\omega^{*}_{j}\rvert^{p-1}]^{2}+p(p-1)\psi^{\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)[\mbox{sgn}(\omega^{*}_{j})]^{2}\lvert\omega^{*}_{j}\rvert^{p-2}$
$\displaystyle=p^{2}\psi^{\prime\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)\lvert\omega^{*}_{j}\rvert^{2p-2}+p(p-1)\psi^{\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)\lvert\omega^{*}_{j}\rvert^{p-2}.$

This completes the proof of the lemma.

We will now show the boundedness of the sequence of Lagrange multiplier
vectors $\\{(\zeta^{k},\eta^{k})\\}$ in the following lemma. The proof follows
the same ideas as those in (OTKW21, , Proposition 9), but is nonetheless
presented to show specifically why the result holds when a general smoothing
function is used.

######  Lemma 3.5.

Suppose that Assumptions (B1)-(B4) and (C1)-(C4) hold. Let
$\\{(\zeta^{k},\eta^{k})\\}\subseteq\Re^{n}\times\Re^{r}$ be a sequence of the
accompanying Lagrange multiplier vectors generated by Algorithm 1. Then
$\\{(\zeta^{k},\eta^{k})\\}$ is bounded.

Proof: For convenience, denote

$\xi^{k}:=((\zeta^{k})^{T},(\eta^{k})^{T})^{T},\quad\hat{\zeta}^{k}:=\frac{\zeta^{k}}{\|\xi^{k}\|},\quad\hat{\eta}^{k}:=\frac{\eta^{k}}{\|\xi^{k}\|}$

for each $k$. We prove by contradiction that the sequence
$\\{(\xi^{k},\eta^{k})\\}$ is bounded. Without loss of generality, we may
assume that

$\|\xi^{k}\|\to\infty,\
\lim_{k\to\infty}\frac{\xi^{k}}{\|\xi^{k}\|}=\hat{\xi}^{*},$

where $\hat{\xi}^{*}:=((\zeta^{*})^{T},(\eta^{*})^{T})^{T}$ with $\zeta^{*}$
and $\eta^{*}$ are accumulation points of $\\{\hat{\zeta}^{k}\\}$ and
$\\{\hat{\eta}^{k}\\}$, respectively. By Assumption (B2), we may suppose
without loss of generality that
$\lim_{k\to\infty}(\omega^{k},\lambda^{k})=(\omega^{*},\lambda^{*})$. Dividing
by $\|\xi^{k}\|$ both sides of (21),(22),(23) and (25) evaluated at
$(\omega,\lambda,\zeta,\eta)=(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})$, we
have for each $k$ the following equations:

$\frac{(\nabla
f(\omega^{k}))_{j}}{\|\xi^{k}\|}+\left(\nabla^{2}_{\omega\omega}G(\omega^{k},\bar{\lambda}^{k})\hat{\zeta}^{k}\right)_{j}+\lambda_{1}^{k}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\hat{\zeta}^{k}_{j}=\frac{(\varepsilon_{1}^{k-1})_{j}}{\|\xi^{k}\|}\
(j=1,2,\dots,n),$ (50)
$\nabla\varphi_{\mu_{k-1}}(\omega^{k})^{T}\hat{\zeta}^{k}-\hat{\eta}^{k}_{1}=\frac{\varepsilon_{2}^{k-1}}{\|\xi^{k}\|},$
(51) $\nabla
R_{j}(\omega^{k})^{T}\hat{\zeta}^{k}-\hat{\eta}^{k}_{j}=\frac{(\varepsilon_{3}^{k-1})_{j}}{\|\xi^{k}\|}\
(j=2,3,\dots,r),$
$\lambda_{j}^{k}\hat{\eta}^{k}_{j}\leq\sum_{i=1}^{r}\lambda_{j}^{k}\hat{\eta}^{k}_{j}=\frac{\varepsilon_{5}^{k-1}}{\|\xi^{k}\|},\quad\lambda_{j}^{k}\geq
0,\quad\hat{\eta}^{k}_{j}\geq 0\quad(j=1,2,3,\dots,r),$ (52)

where
$(\omega,\lambda,\zeta,\eta)=(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})$ and
$(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4},\varepsilon_{5})=(\varepsilon^{k-1}_{1},\varepsilon^{k-1}_{2},\varepsilon^{k-1}_{3},\varepsilon^{k-1}_{4},\varepsilon^{k-1}_{5})$.
Since $\lim_{k\to\infty}\frac{\varepsilon_{l}^{k-1}}{\|\xi^{k}\|}=0$,
$l=1,2,3,5$, and $\lim_{k\to\infty}\lambda^{k}_{1}=\lambda^{*}_{1}>0$ by
Assumption (B1), letting $k\to\infty$ in inequality (52) gives

$\hat{\eta}^{*}_{1}=0,$ (53)

and similarly,

$\hat{\eta}^{*}_{j}=0\ (j\notin
I(\lambda^{*}))~{}~{}\mbox{and}~{}~{}\hat{\eta}^{*}_{j}\geq 0\ (j\in
I(\lambda^{*})).$ (54)

Since $\|\hat{\xi}^{*}\|=1$, we get from (54) that

$1=\|\hat{\zeta}^{*}\|^{2}+\|\hat{\eta}^{*}\|^{2}=\|\hat{\zeta}^{*}\|^{2}+\sum_{j\in
I(\lambda^{*})}\lvert\hat{\eta}^{*}_{j}\rvert^{2}.$ (55)

Meanwhile, since $G$ is twice continuously differentiable and $f$ is
continuously differentiable, then both
$\\{\nabla^{2}_{\omega\omega}G(\omega^{k},\bar{\lambda}^{k})\hat{\zeta}^{k}\\}$
and $\left\\{\frac{\nabla f(\omega^{k})}{\|\xi^{k}\|}\right\\}$ are bounded.
This implies the boundedness of
$\\{\lambda_{1}^{k}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\hat{\zeta}^{k}_{j}\\}$
for each $j$ by (50). Consequently, we obtain
$\lim_{k\to\infty}\hat{\zeta}^{k}_{j}=0$ for $j\in I(\omega^{*})$ by using
Assumption (B1) and Lemma 3.3. That is,

$\hat{\zeta}^{*}_{j}=0\ (j\in I(\omega^{*})).$ (56)

By invoking Assumption (B1) and letting $k\to\infty$ in (24), it is clear that

$\lim_{k\to\infty}\lvert\nabla(\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\rvert=\frac{\lvert
F_{j}(\omega^{*},\bar{\lambda}^{*})\rvert}{\lambda_{1}^{*}},$ (57)

where $F_{j}(\omega^{*},\bar{\lambda}^{*})$ is given by $\eqref{Partial-G}$.
This together with (56) gives us

$\lim_{k\to\infty}\sum_{j\in
I(\omega^{*})}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\hat{\zeta}^{k}_{j}=0.$
(58)

On the other hand, we get

$\lim_{k\to\infty}\nabla\varphi_{\mu_{k-1}}(\omega^{k})^{T}\hat{\zeta}^{k}=0$

by letting $k\to\infty$ in (51) and by using (53). Together with (58), it
yields

$\displaystyle
0=\lim_{k\to\infty}\nabla\varphi_{\mu_{k-1}}(\omega^{k})^{T}\hat{\zeta}^{k}$
$\displaystyle=\lim_{k\to\infty}\left(\sum_{j\in
I(\omega^{*})}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\hat{\zeta}^{k}_{j}+\sum_{j\notin
I(\omega^{*})}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\hat{\zeta}^{k}_{j}\right)$
$\displaystyle=\lim_{k\to\infty}\sum_{j\notin
I(\omega^{*})}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\hat{\zeta}^{k}_{j}$
$\displaystyle=\sum_{j\notin I(\omega^{*})}p{\rm
sgn}(\omega^{*}_{j})\lvert\omega^{*}_{j}\rvert^{p-1}\psi^{\prime}(\lvert\omega^{*}_{j}\rvert^{p})\hat{\zeta}^{*}_{j},$

where the last equality holds by (47). With this equation, we may employ
similar techniques as in the proof of (OTKW21, , Proposition 9) to derive the
equation

$\sum_{j\notin
I(\omega^{*})}\hat{\zeta}^{*}_{j}\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})-\sum_{j\in
I(\lambda^{*})}\hat{\eta}^{*}_{j}\nabla_{(\bar{\omega},\lambda)}\lambda_{j}(\omega^{*},\lambda^{*})=0,$
(59)

where $\bar{\omega}:=(\omega_{j})_{j\notin I(\omega^{*})}$ and
$\Phi_{j}(\omega,\lambda)$ $(j\notin I(\omega^{*}))$ are as defined in
Assumption (B4). On the other hand, by Lemma 3.1, we can find a vector
$d\in\Re^{n-\lvert I(\omega^{*})\rvert+r}$ such that (28) and (29) hold. From
(59), we have

$\sum_{j\notin
I(\omega^{*})}\hat{\zeta}^{*}_{j}\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})^{T}d-\sum_{j\in
I(\lambda^{*})}\hat{\eta}^{*}_{j}\nabla_{(\bar{\omega},\lambda)}\lambda_{j}(\omega^{*},\lambda^{*})^{T}d=0.$

Together with equation (28), we obtain

$\sum_{j\in
I(\lambda^{*})}\hat{\eta}^{*}_{j}\nabla_{(\bar{\omega},\lambda)}\lambda_{j}(\omega^{*},\lambda^{*})^{T}d=0.$

Consequently, we have from (54) and (29) that $\hat{\eta}_{j}^{*}=0$ for all
$j\in I(\lambda^{*})$. In turn, (59) implies that

$\sum_{j\notin
I(\omega^{*})}\hat{\zeta}^{*}_{j}\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})=0.$

Since
$\\{\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})\\}_{j\notin
I(\omega^{*})}$ is linearly independent by Lemma 3.1, then
$\hat{\zeta}_{j}^{*}=0$ for all $j\notin I(\omega^{*})$. Together with (56),
we have $\hat{\zeta}^{*}=0$ which in turn implies that
$\|\hat{\zeta}^{*}\|^{2}+\sum_{j\in
I(\lambda^{*})}\lvert\hat{\eta}^{*}_{j}\rvert=0$. This contradicts (55).
Therefore, the sequence $\\{(\zeta^{k},\eta^{k})\\}$ is bounded.

By the above lemma, we know that accumulation points of
$\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}$ exist. We now show that
such point satisfies condition (9)-(10).

######  Lemma 3.6.

Suppose that Assumptions (B1)-(B4) and (C1)-(C4) hold. Then, an accumulation
point of the sequence $\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}$
generated by Algorithm 1 satisfies

  * (i)

$\zeta_{j}^{*}=0$ for all $j\in I(\omega^{*})$ satisfies condition (10).

  * (ii)

$p\sum_{j\notin I(\omega^{*})}{\rm
sgn}(\omega^{*}_{j})\lvert\omega^{*}_{j}\rvert^{p-1}\psi^{\prime}(\lvert\omega^{*}_{j}\rvert^{p})\zeta^{*}_{j}=\eta^{*}_{1}$
satisfies condition (9).

Proof: The proof essentially follows ideas similar to Lemma 3.5. Indeed, from
Lemma 3.5 and Assumption (B2), the sequence
$\\{(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})\\}$ is bounded. Hence, without
loss of generality, we assume that

$\lim_{k\to\infty}(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})=(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*}).$
(60)

Since $G$ is twice continuously differentiable and $f$ is continuously
differentiable, then
$\\{\lambda_{1}^{k}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}\zeta^{k}_{j}\\}$
is bounded for each $j$ by using (60) and equation (21). From Lemma 3.3 and
Assumption (B1), we obtain that $\lim_{k\to\infty}\zeta^{k}_{j}=0$ for $j\in
I(\omega^{*})$, that is, $\zeta^{*}_{j}=0$ for $j\in I(\omega^{*})$. This
proves (i).

We now prove (ii). From (57), we know that
$\\{(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\\}$ is bounded. Hence, we
obtain from item (i) that
$\lim_{k\to\infty}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\zeta^{k}_{j}=0$,
$(j\in I(\omega^{*}))$ and therefore,

$\lim_{k\to\infty}\sum_{j\in
I(\omega^{*})}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\zeta^{k}_{j}=0.$

Together with (47), it follows that

$\displaystyle\lim_{k\to\infty}\nabla\varphi_{\mu_{k-1}}(\omega^{k})^{T}\zeta^{k}$
$\displaystyle=\lim_{k\to\infty}\sum_{j\notin
I(\omega^{*})}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}\zeta^{k}_{j}$
$\displaystyle=\sum_{j\notin I(\omega^{*})}p{\rm
sgn}(\omega^{*}_{j})\lvert\omega^{*}_{j}\rvert^{p-1}\psi^{\prime}(\lvert\omega^{*}_{j}\rvert^{p})\zeta^{*}_{j}.$

Using this and taking the limit in (22), we obtain item (ii). The proof is
complete.

We now establish our last lemma which will later be used to prove that
accumulation points of the sequence generated by the algorithm satisfy
conditions (7)-(8). The following result, similar to Lemma 3.4, is in fact a
general formula. That is, the sequence $\\{(\omega^{k},\mu_{k-1})\\}$ need not
be generated by Algorithm 1.

######  Lemma 3.7.

Suppose that Assumption (C1)-(C2) hold. Let
$\\{(\omega^{k},\mu_{k})\\}\subseteq\Re^{n}\times\Re_{++}$ be an arbitrary
sequence converging to $(\omega^{*},0)$. Then, we have

  * (i)

$\lim_{k\to\infty}W_{k}\nabla\varphi_{\mu_{k-1}}(\omega^{k})=p\lvert
W^{*}\rvert^{p}\psi^{\prime}(\lvert\omega^{*}\rvert^{p})$.

  * (ii)

$\lim_{k\to\infty}W^{2}_{k}\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k})=p(p-1){\rm
diag}(\lvert
W^{*}\rvert^{p}\psi^{\prime}(\lvert\omega^{*}\rvert^{p}))+p^{2}{\rm
diag}(\lvert W^{*}\rvert^{2p}\psi^{\prime\prime}(\lvert\omega^{*}\rvert^{p}))$

where $W_{k}:={\rm diag}(\omega^{k})$ for each $k$, $\lvert
W^{*}\rvert^{p}:={\rm diag}(\lvert\omega^{*}\rvert^{p})$, and $\lvert
W^{*}\rvert^{2p}:={\rm diag}(\lvert\omega^{*}\rvert^{2p})$.

Proof: We show first that (i) holds. To this end, we need to show that for
each $j\in\\{1,\dots,n\\}$, we have

$\displaystyle\lim_{k\to\infty}\omega^{k}_{j}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}$
$\displaystyle=\lim_{k\to\infty}p\omega_{j}^{k}\psi^{\prime}\left([\phi(\mu_{k-1},\omega_{j}^{k})]^{p}\right)\phi^{\prime}(\mu_{k-1},\omega_{j}^{k})[\phi(\mu_{k-1},\omega_{j}^{k})]^{p-1}{}$
$\displaystyle=p\lvert\omega^{*}_{j}\rvert^{p}\psi^{\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right),$
(61)

where the first equality holds because of (34). In the case that $j\notin
I(\omega^{*})$, the above can be easily verified by using equation (47). Now,
suppose that $j\in I(\omega^{*})$, i.e. $\omega_{j}^{*}=0$. Let
$K=\\{1,2,\dots\\}$ and as in Lemma 3.3, define

$U^{j}_{1}:=\\{k\in K~{}\lvert~{}\omega_{j}^{k}=0\\},\quad\text{and}\quad
U^{j}_{2}:=\\{k\in K~{}\lvert~{}\omega_{j}^{k}\neq 0\\}.$

It is clear that

$\lim_{k\in
U_{1}^{j}\to\infty}\omega^{k}_{j}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}=0.$
(62)

For $k\in U^{j}_{2}$, we have from Proposition 2.7(b) that
$\phi(\mu_{k-1},\omega^{k}_{j})\geq\lvert\omega^{k}_{j}\rvert$, which gives us

$[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}\leq\lvert\omega^{k}_{j}\rvert^{p-1}$

for any $p\in(0,1]$. Together with Proposition 2.7(c), we have

$\lvert\omega^{k}_{j}\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}\rvert\leq\lvert\omega^{k}_{j}\rvert.\lvert\omega^{k}_{j}\rvert^{p-1}=\lvert\omega^{k}_{j}\rvert^{p}\to\lvert\omega^{*}_{j}\rvert^{p}=0\
\mbox{as}\ k\in U_{2}^{j}\to\infty.$

Since $0<\psi^{\prime}(t)\leq\alpha$ by Assumption (A), we get

$\lim_{k\in
U_{2}^{j}\to\infty}\omega^{k}_{j}(\nabla\varphi_{\mu_{k-1}}(\omega^{k}))_{j}=0.$
(63)

Equations (62) and (63) confirm the desired equation (61).

Next, we show that (ii) holds. It is enough to show that for each
$j\in\\{1,\dots,n\\}$, we have

$\lim_{k\to\infty}(\omega^{k}_{j})^{2}(\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))_{jj}=p^{2}\psi^{\prime\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)\lvert\omega^{*}_{j}\rvert^{2p}+p(p-1)\psi^{\prime}\left(\lvert\omega^{*}_{j}\rvert^{p}\right)\lvert\omega^{*}_{j}\rvert^{p}.$
(64)

The case when $j\notin I(\omega^{*})$ easily follows from (48). Now, suppose
that $j\in I(\omega^{*})$. Since $\phi$ is twice continuously differentiable,
we may consider the Maclaurin expansion of $\phi$:

$\displaystyle\phi(\mu_{k-1},\omega^{k}_{j})=\phi(\mu_{k-1},\lvert\omega^{k}_{j}\rvert)$
$\displaystyle=\phi(\mu_{k-1},0)+\phi^{\prime}(\mu_{k-1},0)\lvert\omega^{k}_{j}\rvert+\frac{\phi^{\prime\prime}(\mu_{k-1},\gamma^{k}_{j})}{2}\lvert\omega^{k}_{j}\rvert^{2}$
$\displaystyle=\phi(\mu_{k-1},0)+\frac{\phi^{\prime\prime}(\mu_{k-1},\gamma^{k}_{j})}{2}\lvert\omega^{k}_{j}\rvert^{2}$

where $\gamma^{k}_{j}\in[0,\lvert\omega^{k}_{j}\rvert)$, the first equality
holds by Assumption (C1) and (17), and the last equality follows from (30).
Now, since $\phi(\mu_{k-1},0)>0$,

$\displaystyle\phi^{\prime\prime}(\mu_{k-1},\gamma^{k}_{j})\lvert\omega^{k}_{j}\rvert^{2}$
$\displaystyle=2(\phi(\mu_{k-1},\omega^{k}_{j})-\phi(\mu_{k-1},0))\leq
2\phi(\mu_{k-1},\omega^{k}_{j}).$

Then

$\displaystyle[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}\phi^{\prime\prime}(\mu_{k-1},\gamma^{k}_{j})(\omega^{k}_{j})^{2}$
$\displaystyle\leq 2[\phi(\mu_{k-1},\omega^{k}_{j})]^{p}\to
2\lvert\omega^{*}_{j}\rvert^{p}=0\ \mbox{as}\ k\to\infty.$ (65)

On the other hand, since $\rho$ is nonincreasing on $[0,\infty)$,

$\phi^{\prime\prime}(\mu_{k-1},\gamma^{k}_{j})=\frac{2}{\mu_{k-1}}\rho\left(\frac{\gamma^{k}_{j}}{\mu_{k-1}}\right)\geq\frac{2}{\mu_{k-1}}\rho\left(\frac{\lvert\omega^{k}_{j}\rvert}{\mu_{k-1}}\right)=\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j}).$
(66)

It follows from (65) and (66) that

$\lim_{k\to\infty}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}\phi^{\prime\prime}(\mu_{k-1},\omega^{k}_{j})(\omega^{k}_{j})^{2}=0.$
(67)

Next, we will show that

$\lim_{k\to\infty}(\omega_{j}^{k})^{2}[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-2}=0.$
(68)

Using again Proposition 2.7(b), we get
$[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-2}\leq\lvert\omega^{k}_{j}\rvert^{p-2}$.
Then

$(\omega_{j}^{k})^{2}[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})]^{2}[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-2}\leq\lvert\omega^{k}_{j}\rvert^{p}.$

Letting $k\to\infty$, we obtain the desired limit (68). As a consequence of
(68), we obtain

$\lim_{k\to\infty}(\omega_{j}^{k})^{2}[\phi^{\prime}(\mu_{k-1},\omega^{k}_{j})[\phi(\mu_{k-1},\omega^{k}_{j})]^{p-1}]^{2}=0.$
(69)

Equation (64) now follows by invoking Assumption (A), using the formula (35)
and the derived limits (67), (68) and (69). This completes the proof of (ii).

Having derived all the necessary lemmas, we can now prove our main result.

Proof: [Proof of Theorem 3.1] Let
$(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*})$ be an arbitrary accumulation
point of the sequence $(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})$. By taking
a subsequence, if necessary, we may assume without loss of generality that

$\lim_{k\to\infty}(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})=(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*}).$

We need to show that $(\omega^{*},\lambda^{*},\zeta^{*},\eta^{*})$ satisfies
SB-KKT conditions (7)-(12). Indeed, multiplying by $W^{2}_{k}$ and $W_{k}$
both sides of (21) and (24) evaluated at
$(\omega,\lambda,\zeta,\eta)=(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})$,
respectively, we have

$W^{2}_{k}\nabla
f(\omega^{k})+W^{2}_{k}(\nabla^{2}_{\omega\omega}G(\omega^{k},\bar{\lambda}^{k})+\lambda_{1}^{k}\nabla^{2}\varphi_{\mu_{k-1}}(\omega^{k}))\zeta^{k}=W^{2}_{k}\varepsilon^{k-1}_{1},$
$W_{k}\nabla_{\omega^{k}}G(\omega^{k},\bar{\lambda}^{k})+\lambda_{1}^{k}W_{k}\nabla\varphi_{\mu_{k-1}}(\omega^{k})=W_{k}\varepsilon^{k-1}_{4}$

where $W_{k}=\mbox{diag}(\omega^{k})$,
$(\omega,\lambda,\zeta,\eta)=(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})$ and

$(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4},\varepsilon_{5})=(\varepsilon^{k-1}_{1},\varepsilon^{k-1}_{2},\varepsilon^{k-1}_{3},\varepsilon^{k-1}_{4},\varepsilon^{k-1}_{5}).$

Recall that $G$ is twice continuously differentiable and $f$ is continuously
differentiable. Thus, letting $k\to\infty$ in the above two equations and
using Lemma 3.7, we obtain that

$W^{2}_{*}\nabla f(\omega^{*})+H(\omega^{*},\lambda^{*})\zeta^{*}=0,$
$W_{*}\nabla_{\omega^{*}}G(\omega^{*},\bar{\lambda}^{*})+\lambda_{1}^{*}W_{*}\nabla\varphi_{\mu_{k-1}}(\omega^{*})=0,$

where

$H(\omega,\lambda)=W^{2}\nabla^{2}_{\omega}G(\omega,\bar{\lambda})+\lambda_{1}p(p-1){\rm
diag}(\lvert
W\rvert^{p}\psi^{\prime}(\lvert\omega\rvert^{p}))+\lambda_{1}p^{2}{\rm
diag}(\lvert W\rvert^{2p}\psi^{\prime\prime}(\lvert\omega\rvert^{p})).$

The conditions (9) and (10) can be derived by Lemma 3.6. Let $k\to\infty$ in
(23) and (25) with
$(\omega,\lambda,\zeta,\eta)=(\omega^{k},\lambda^{k},\zeta^{k},\eta^{k})$, we
obtain the conditions (11) and (12). This completes the proof of Theorem 3.1.

##  4 Numerical results

###  4.1 Examples of Smoothing Functions

A key aspect in successful numerical implementations of a smoothing algorithm
is the choice of the approximating functions. Here, we enumerate six smoothing
functions that we will use in our numerical simulations.

There are many density functions commonly used and called kernel functions in
statistics (see also CM96 ). Some density functions satisfying (13) are given
as follows.

$\displaystyle\rho_{1}(x)$ $\displaystyle:=$
$\displaystyle\left\\{\begin{array}[]{clc}\frac{35}{32}(1-x^{2})^{3}&{\rm
if}\;\lvert x\rvert\leq 1,\\\ 0&\mbox{otherwise}.\end{array}\right.$
$\displaystyle\rho_{2}(x)$ $\displaystyle:=$
$\displaystyle\left\\{\begin{array}[]{clc}\frac{15}{16}(1-x^{2})^{2}&\mbox{if}\;\lvert
x\rvert\leq 1,\\\ 0&\mbox{otherwise}.\end{array}\right.$
$\displaystyle\rho_{3}(x)$ $\displaystyle:=$
$\displaystyle\left\\{\begin{array}[]{clc}\frac{3}{4}(1-x^{2})&\mbox{if}\;\lvert
x\rvert\leq 1,\\\ 0&\mbox{otherwise},\end{array}\right.$
$\displaystyle\rho_{4}(x)$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}\quad\forall x\in\Re.$
$\displaystyle\rho_{5}(x)$ $\displaystyle:=$
$\displaystyle\frac{e^{-x}}{(1+e^{-x})^{2}}.$ $\displaystyle\rho_{6}(x)$
$\displaystyle:=$ $\displaystyle\frac{1}{(x^{2}+4)^{\frac{3}{2}}}.$

Following the discussion in Section 2.3, the corresponding smoothing functions
of $\lvert x\rvert$ are given as follows:

$\displaystyle\phi_{1}(\mu,x)$ $\displaystyle:=$
$\displaystyle\left\\{\begin{array}[]{clc}-\frac{5x^{8}}{128\mu^{7}}+\frac{7x^{6}}{32\mu^{5}}-\frac{35x^{4}}{64\mu^{3}}+\frac{35x^{2}}{32\mu}+\frac{35\mu}{128}&\mbox{if}\;\lvert
x\rvert\leq\mu,\\\ \lvert x\rvert&\mbox{if}\;\lvert
x\rvert>\mu.\end{array}\right.$ $\displaystyle\phi_{2}(\mu,x)$
$\displaystyle:=$
$\displaystyle\left\\{\begin{array}[]{clc}\frac{x^{6}}{16\mu^{5}}-\frac{5x^{4}}{16\mu^{3}}+\frac{15x^{2}}{16\mu}+\frac{5\mu}{16}&\mbox{if}\;\lvert
x\rvert\leq\mu,\\\ \lvert x\rvert&\mbox{if}\;\lvert
x\rvert>\mu.\end{array}\right.$ $\displaystyle\phi_{3}(\mu,x)$
$\displaystyle:=$
$\displaystyle\left\\{\begin{array}[]{clc}-\frac{x^{4}}{8\mu^{3}}+\frac{3x^{2}}{4\mu}+\frac{3\mu}{8}&\mbox{if}\;\lvert
x\rvert\leq\mu,\\\ \lvert x\rvert&\mbox{if}\;\lvert
x\rvert>\mu.\end{array}\right.$ $\displaystyle\phi_{4}(\mu,x)$
$\displaystyle:=$ $\displaystyle x{\rm
erf}\left(\frac{x}{\sqrt{2}\mu}\right)+\sqrt{\frac{2}{\pi}}\mu
e^{-\frac{x^{2}}{2\mu^{2}}},$ $\displaystyle\phi_{5}(\mu,x)$ $\displaystyle:=$
$\displaystyle\mu\begin{bmatrix}\log\begin{pmatrix}1+e^{-\frac{x}{\mu}}\\\
\end{pmatrix}+\log\begin{pmatrix}1+e^{\frac{x}{\mu}}\\\ \end{pmatrix}\\\
\end{bmatrix}.$ $\displaystyle\phi_{6}(\mu,x)$ $\displaystyle:=$
$\displaystyle\sqrt{4\mu^{2}+x^{2}}.$

Here, the error function is defined by

$\mbox{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-u^{2}}du\quad\forall
x\in\Re.$

The graphs of $\lvert x\rvert$ and $\phi_{i}(\mu,x)$, $i=1,2,\dots,6$ with
$\mu=0.25$ are illustrated in Figure 1. From the graphs, we infer the
following inequality relating the smoothing functions:

$\lvert
x\rvert\leq\phi_{1}(\mu,x)\leq\phi_{2}(\mu,x)\leq\phi_{3}(\mu,x)\leq\phi_{4}(\mu,x)\leq\phi_{5}(\mu,x)\leq\phi_{6}(\mu,x).$

It is not difficult to show that the relation $\lvert
x\rvert\leq\phi_{1}(\mu,x)\leq\phi_{2}(\mu,x)\leq\phi_{3}(\mu,x)$, while the
proof of the relation
$\phi_{3}(\mu,x)\leq\phi_{4}(\mu,x)\leq\phi_{5}(\mu,x)\leq\phi_{6}(\mu,x)$ can
be found in SNC19 . On the other hand, the graphs of the corresponding
smoothing functions for $\lvert x\rvert^{p}$ where $p\in(0,1]$ is shown in
Figure 2.

Figure 1: Graph of $\lvert x\rvert$ and $\phi_{i}(\mu,x)$, $i=1,2,\dots,6$
with $\mu=0.25$.
Figure 2: Graph of $\lvert x\rvert^{p}$ and $(\phi(\mu,x))^{p}$,
$i=1,2,\dots,6$ with $\mu=0.25$ and $p=0.5$.

The smooth approximation $\phi(\mu,x):=\sqrt{\mu^{2}+x^{2}}$ is the function
used in OTKW21  for their smoothing algorithm for (3) with
$R_{1}(\omega):=\sum_{i=1}^{n}\lvert\omega_{i}\rvert^{p}$ $(0<p\leq 1)$. This
function is simply $\phi_{6}(\frac{\mu}{2},x)$ and can in fact be obtained by
choosing the density function $\rho(t)=\frac{1}{2}(x^{2}+1)^{-3/2}$.

In this paper, we consider the six functions above and determine which
approximation is the best suitable in solving (3) with $R_{1}(\omega)$
satisfying Assumption (A). It is easy to check that the six density functions
as above satisfy Assumptions (C1)-(C2). Condition (C3), on the other hand,
holds by choosing $c=4$, and $r=2$. Indeed,

$1-\frac{4}{4+S^{2}}\leq\sqrt{1-\frac{4}{4+S^{2}}}=2\int\limits_{0}^{S}\rho_{3}(s)~{}ds\leq
2\int\limits_{0}^{S}\rho_{i}(s)~{}ds~{}~{}\forall i=1,\dots,6.$

According to Assumption (C4), only the functions $\rho_{4}$, $\rho_{5}$ and
$\rho_{6}$ can be used (theoretically) for the case $p=1$.

###  4.2 Simulations

We compare efficiency of the smoothing functions $\phi_{i}\ (i=1,2,\ldots,6)$
by means of numerical simulation. The program is coded in MATLAB R2019a and
run on a machine with Intel(R) Core(TM) i7-8565U<EMAIL_ADDRESS>and 16.0 GB RAM.

####  4.2.1 Problem with an Elastic-Net-type regularizer

We solve the following bilevel problem arising from squared linear regression
using an Elastic-Net-type regularizer:

$\displaystyle\begin{array}[]{rll}\displaystyle{\min_{\omega,\lambda}}&&\displaystyle{\frac{1}{2}\|A_{1}{\omega}-b_{1}\|_{2}^{2}}\\\
\mbox{s.t. }&&\omega\in\displaystyle{\mathop{\rm
argmin}_{\hat{\omega}}}\left\\{\frac{1}{2}\|A_{2}\hat{\omega}-b_{2}\|_{2}^{2}+{\lambda_{1}\|\hat{\omega}\|_{p}^{p}}+\lambda_{2}\|\hat{\omega}\|_{2}^{2}\right\\}\\\
&&(\lambda_{1},\lambda_{2})\geq 0,\end{array}$ (79)

where $A_{i}\in\Re^{m_{i}\times n}$ and $b_{i}\in\Re^{m_{i}}$ for
$i\in\\{1,2\\}$. We produce 20 synthetic problems for each case of
$(n,m_{1},m_{2})=(500,300,300),(500,1000,1000)$ by

$\displaystyle A_{i}:=\texttt{rand}(m_{i},n),\ \begin{bmatrix}b_{1}\\\
b_{2}\end{bmatrix}:=\begin{bmatrix}A_{1}*\theta\\\
A_{2}*\theta+0.01*\left(2*\texttt{rand}\left(m_{2},1\right)-\texttt{ones}\left(m_{2},1\right)\right)\end{bmatrix},$
$\displaystyle\theta:=10*\left[\texttt{zeros}(\lfloor\frac{n}{2}\rfloor,1)\mathchar
24635\relax\;\texttt{rand}(n-\lfloor\frac{n}{2}\rfloor,1)\right]$

with rand, ones, and zeros being MATLAB commands, and apply Algorithms 1 with
the smoothing functions $\phi_{i}\ (i=1,2,\ldots,6)$ to these problems. The
random number generator is initialized by $\texttt{rng}(0.111111)$. Moreover,
the parameter $p$ in the $\ell_{p}$-regularizer is set to $p=1,0.5$ for each
problem.

In order to compute a KKT point of the smoothed subproblem for (79) in Step 1
of Algorithm 1, we utilize the MATLAB solver fmincon with
“MaxIterations$=10^{7}$” and opt for the interior-point method as an algorithm
that runs within fmincon.222 One may think of using the SQP method in place of
the interior-point method. In the preliminary numerical experiments, fmincon
using the interior point method (fminconIPM) performed better than the one
with the SQP method (fminconSQP). Indeed, we observed that fminconIPM solved
each smoothed subproblem for (79) more stably than fminconSQP. We also
observed that fminconIPM was more likely to find a positive solution
$\lambda_{1}$, leading to a sparse solution. This result may be peculiar to an
interior point method that reaches a solution from the interior of a feasible
region.  As a starting point of fmincon, we choose
$\omega=100*\texttt{ones}(n,1)$ at the first iteration $k=0$ and then use the
previous iteration point $(\omega^{k-1},\lambda^{k-1})$ after the second
iteration, i.e., for $k\geq 1$. We set $(\mu_{0},\beta_{1})=(1,0.8)$ in Step 0
for Algorithm 1. To compute the value of the function ${\rm erf}$ in
$\phi_{4}$, we utilize the MATLAB function erf.

In light of the SB-KKT conditions (7)-(12) and the value of the smoothing
parameter $\mu$, we terminate the algorithm when either one of the following
criteria is met:

  1. 1.

The norms of the residuals of the equations in (7)-(12) are smaller than
$10^{-2}$. As for condition (12), we observe only the value of
$(\lambda^{\ast})^{\top}\eta^{\ast}$ because the nonnegativity of
$(\lambda^{\ast},\eta^{\ast})$ is assured in the interior-point method. To
estimate the index set $I(\omega^{\ast})$ in conditions (9) and (10), we
regard $\omega^{k}_{i}$ as zero if $\lvert\omega^{k}_{i}\rvert\leq 10^{-5}$.

  2. 2.

$\mu_{k+1}\leq 10^{-8}$.

The obtained results are summarized in Tables 1-4, in which each column is
described as follows. Here, the averages are taken over the set of problems
that are counted in success($\%$).

$i$: | the smoothing function $\phi_{i}$
---|---
obj: | averaged objective value at the resulting solution
sbkkt: | averaged residual of the SB-KKT conditions
$\mu_{\rm end}$: | averaged value of $\mu_{k+1}$ on termination
${\rm sparsity}(\%)$: | averaged ratio of zero elements of the resulting solution $\omega$,
| in which each element $w_{i}$ is counted as zero if
$\lvert\omega_{i}\rvert\leq 10^{-5}$
ave.time(s): | averaged spent time in seconds
ite: | averaged number of iterations
success($\%$): | ratio of problems for which SB-KKT points are computed successfully
| in the sense that the first termination condition in the above is satisfied

The best values for obj, sbkkt, ${\rm sparsity}(\%)$, and ave.time(s) are
displayed in bold in the tables. In addition, the smallest ones for $\mu_{\rm
end}$ and ite are also in bold. The obtained results with $\phi_{5}$ were
excluded from the tables because the overflow often occurred when computing
its gradient as $\mu$ gets smaller and it did not work well. Now, the
following insights are obtained from the numerical results:

##### Comparison of results with $p=0.5,1$:

From ${\rm sparsity(\%)}$, we see that $\ell_{0.5}$ tends to attain sparser
solutions than $\ell_{1}$. For example, by Tables 3 and 4, ${\rm
sparsity(\%)}$ for $p=1$ takes values less than 10%, whereas they are more
than 60% for $p=0.5$. This is by virtue of the nonconvexity of $\ell_{p}$ with
$p<1$. Moreover, $\ell_{0.5}$ tends to attain solutions with better objective
values than $\ell_{1}$. Indeed, values of obj in Table 3 with $p=1$ are around
$3.5\cdot 10^{-3}$, while those in Table 4 with $p=0.5$ are around $7.2\cdot
10^{-4}$. Similar observations can be gathered from Tables 1 and 2.

##### Comparison of results with $m_{2}=300,1000$:

The problems with $m_{2}<n$ can be related to the problem of finding a sparse
solution of underdetermined linear system. Such kind of problems are often
regarded more intractable than those with $m_{2}\geq n$, as illustrated by the
obtained numerical results: Any ave.time(s) are more than 450 seconds in
Tables 1 and 2 with $m_{2}=300<n=500$, whereas they are less than 400 seconds
in Tables 3 and 4 with $m_{2}=1000>n$.

##### Comparison of results of the four smoothing functions:

In view of the objective values and residual of the SB-KKT conditions, the
qualities of the resulting solutions are comparable. Meanwhile, substantial
differences are observed in tendency of the smoothing parameter $\mu$ on
termination, i.e., $\mu_{end}$. The function $\phi_{6}$ takes the smallest
$\mu_{end}$ throughout the tables, while $\phi_{1}$ the largest. A possible
reason for this situation is that as illustrated by Figures 1 and 2,
$\phi_{6}(x)$ deviates from $\lvert x\rvert$ most, while $\phi_{1}(x)$ fits to
it the best.

On the other hand, from ave.time(s), we observe that $\phi_{4}$ attains the
smallest cpu time. In particular, cpu-time per a single iteration for
$\phi_{4}$ also tends to be the smallest. Indeed, for example, from Table 4,
$\frac{ave.time(s)}{ite}$ is around 7.5 seconds for $\phi_{4}$ while more than
8 seconds for the others. This fact may suggest that smoothed subproblems
organized with $\phi_{4}$ are easier for fmincon to solve than those with the
other smoothing functions.

Table 1: Averaged results for $(n,m_{2},m_{1},p)=(500,300,300,1)$ $i$ | obj | sbkkt | $\mu_{\rm end}$ | ${\rm sparsity}(\%)$ | ave.time(s) | ite | success(%)
---|---|---|---|---|---|---|---
1 | 3.77e-03 | 8.06e-03 | 4.07e-04 | 5.17 | 509.95 | 35.2 | 100
2 | 3.75e-03 | 7.75e-03 | 3.65e-04 | 5.17 | 514.40 | 35.7 | 95
3 | 3.74e-03 | 7.19e-03 | 2.97e-04 | 5.08 | 505.11 | 36.6 | 95
4 | 3.76e-03 | 8.16e-03 | 1.44e-04 | 5.31 | 418.04 | 39.8 | 95
6 | 3.77e-03 | 8.10e-03 | 5.55e-05 | 7.96 | 473.83 | 44.1 | 100
Table 2: Averaged results for $(n,m_{2},m_{1},p)=(500,300,300,0.5)$ (Note: Although obj with $i=3$ is remarkably higher than the others in the table, this is caused by one outlier instance. After precluding that outlier, we confirmed that the retaken obj was comparative to the others.) $i$ | obj | sbkkt | $\mu_{\rm end}$ | ${\rm sparsity}(\%)$ | ave.time(s) | ite | success(%)
---|---|---|---|---|---|---|---
1 | 7.88e-04 | 5.56e-03 | 5.78e-04 | 54.29 | 473.49 | 34.6 | 90
2 | 8.32e-04 | 6.16e-03 | 5.17e-04 | 50.68 | 450.14 | 35.3 | 95
3 | 6.64e+02 | 5.32e-03 | 3.44e-04 | 56.60 | 524.97 | 37.8 | 95
4 | 8.34e-04 | 6.21e-03 | 2.06e-04 | 52.61 | 389.90 | 39.3 | 100
6 | 8.38e-04 | 5.94e-03 | 1.31e-04 | 52.07 | 479.38 | 41.7 | 90
Table 3: Averaged results for $(n,m_{2},m_{1},p)=(500,1000,1000,1)$ $i$ | obj | sbkkt | $\mu_{\rm end}$ | ${\rm sparsity}(\%)$ | ave.time(s) | ite | success(%)
---|---|---|---|---|---|---|---
1 | 2.76e-03 | 8.17e-03 | 4.31e-04 | 5.00 | 158.12 | 34.8 | 100
2 | 2.77e-03 | 8.33e-03 | 3.83e-04 | 4.83 | 157.58 | 35.3 | 100
3 | 2.79e-03 | 8.73e-03 | 3.27e-04 | 4.71 | 170.79 | 36 | 100
4 | 2.78e-03 | 8.30e-03 | 1.52e-04 | 5.06 | 137.85 | 39.5 | 100
6 | 2.81e-03 | 8.50e-03 | 7.27e-05 | 6.49 | 154.41 | 42.8 | 100
Table 4: Averaged results for $(n,m_{2},m_{1},p)=(500,1000,1000,0.5)$ $i$ | obj | sbkkt | $\mu_{\rm end}$ | ${\rm sparsity}(\%)$ | ave.time(s) | ite | success(%)
---|---|---|---|---|---|---|---
1 | 7.28e-04 | 6.00e-03 | 4.51e-04 | 60.97 | 338.62 | 35.1 | 100
2 | 7.28e-04 | 6.70e-03 | 3.92e-04 | 61.09 | 325.36 | 35.7 | 100
3 | 7.32e-04 | 6.46e-03 | 3.22e-04 | 60.75 | 329.30 | 36.5 | 95
4 | 7.28e-04 | 6.70e-03 | 1.62e-04 | 60.92 | 298.07 | 39.6 | 100
6 | 7.27e-04 | 6.43e-03 | 8.53e-05 | 61.95 | 343.90 | 42.5 | 100

####  4.2.2 Problems with other regularizers

In this section, we solve problem (79) with the regularizers
$\psi_{2}(\|\omega\|_{p}^{p})$ and $\psi_{3}(\|\omega\|_{p}^{p})$ in place of
$\|\omega\|_{p}^{p}$, where $\psi_{2}$ and $\psi_{3}$ are defined in Appendix
A. We set $a=1$ in $\psi_{2}$ and $\psi_{3}$, and also set $\beta_{1}=0.5$
together with $(n,m_{\rm tr},m_{\rm val})=(500,1000,1000)$. The experimental
settings are the same as in the previous section, and in particular the 20
synthetic problem-data of $A_{i},b_{i}\ (i=1,2)$ are also identical to the
ones used there. The obtained results are summarized in Tables 5 and 6. In
both the cases, more than 95% instances are successfully solved.

In both Tables 5 and 6, we observe that $\mu_{\rm end}$ for $\phi_{6}$ are the
smallest among the five smoothing functions, respectively, as well as the
results in Table 4 with the $\ell_{0.5}$ regularizer. Moreover, those of
$\phi_{1}$ are almost the largest although the exactly largest one in Table 5
is attained by $\phi_{2}$. However, unlike the $\ell_{0.5}$ regularizer,
superiority of $\phi_{4}$ in ave.time(s) is not observed. In comparison of
Tables 4, 5, and 6, the $\ell_{0.5}$ regularizer tends to gain sparser
solutions according to ${\rm sparsity}(\%)$. On the other hand, from ave-
time(s), we see that the spent time for $\psi_{3}$ tends to be the smaller
than those for the $\ell_{0.5}$-regularizer and $\psi_{2}$.

Table 5: Averaged results for $(n,m_{\rm tr},m_{\rm val},p)=(500,1000,1000,0.5)$ using $\psi_{2}$ $i$ | obj | sbkkt | $\mu_{\rm end}$ | ${\rm sparsity}(\%)$ | ave.time(s) | ite | success(%)
---|---|---|---|---|---|---|---
1 | 5.53e-03 | 5.42e-03 | 4.49e-04 | 26.95 | 308.42 | 36.7 | 95
2 | 5.57e-03 | 4.37e-03 | 5.41e-04 | 27.40 | 298.43 | 37.0 | 100
3 | 5.57e-03 | 4.80e-03 | 4.36e-04 | 25.48 | 290.19 | 37.8 | 100
4 | 5.53e-03 | 4.55e-03 | 1.85e-04 | 28.40 | 304.39 | 41.4 | 95
6 | 5.57e-03 | 5.40e-03 | 9.13e-05 | 27.10 | 322.90 | 44 | 100
Table 6: Averaged results for $(n,m_{\rm tr},m_{\rm val},p)=(500,1000,1000,0.5)$ using $\psi_{3}$ $i$ | obj | sbkkt | $\mu_{\rm end}$ | ${\rm sparsity}(\%)$ | ave.time(s) | ite | success(%)
---|---|---|---|---|---|---|---
1 | 5.58e-03 | 4.38e-03 | 5.98e-04 | 26.15 | 285.68 | 36.5 | 100
2 | 5.55e-03 | 4.93e-03 | 3.60e-04 | 28.14 | 295.70 | 37.9 | 100
3 | 5.99e-03 | 5.23e-03 | 3.90e-04 | 25.88 | 304.30 | 38.5 | 100
4 | 5.56e-03 | 5.12e-03 | 1.66e-04 | 27.64 | 321.49 | 41.5 | 100
6 | 5.59e-03 | 4.41e-03 | 1.21e-04 | 27.12 | 313.94 | 43.8 | 100

##  5 Conclusion

This paper considers a class of nonsmooth, possibly nonconvex and non-
Lipschitz regularizers for the best hyperparameter selection problem using a
bilevel programming strategy. The class of regularizers we consider subsumes
the traditional $\ell_{p}$ regularizer, which is the focus of the earlier work
OTKW21 . Moreover, the convergence analysis of the smoothing algorithm
presented in this paper is unified, in the sense that it is not limited to the
chosen smoothing function, unlike the previous work OTKW21  where the analysis
is centered on the selected smoothing function. Finally, we proved our main
convergence result under a milder constraint qualification. More precisely, we
only assumed the Mangasarian-Fromovitz constraint qualification (MFCQ) for our
convergence analysis, which is weaker than the linearly independent constraint
qualification (LICQ) used in OTKW21 . For our numerical simulations, we
compared the performance of six smoothing functions in solving the bilevel
programming problem using different regularizers. Theoretically, we can use
these smoothing functions for all the regularizers considered as their
corresponding density functions satisfy Assumption (C). On the other hand, our
practical experience revealed that the choice of smoothing function is
dependent as well on the chosen regularizer. When using the traditional
$\ell_{p}$ $(p=1,0.5)$ regularizer, the function $\phi_{4}$ has the best
performance in terms of CPU time spent. For the regularizers $\psi_{2}$ and
$\psi_{3}$, the smoothing functions $\phi_{3}$ and $\phi_{1}$, respectively,
have the most competitive performance in terms of CPU time. However, we see
that the $\ell_{p}$ regularizer tends to obtain sparser solutions than the
other two.

## Acknowledgement

This work was conducted while C. T. Nguyen was a research assistant at
National Taiwan Normal University and J. H. Alcantara was a postdoctoral
fellow at National Taiwan Normal University. The work of J.-S. Chen was
supported by the Ministry of Science and Technology, Taiwan.

##  Appendix A Penalty functions that satisfy Assumption (A)

We consider four penalty functions as follows:

$\psi_{1}(t)=t,~{}~{}\psi_{2}(t)=\log(1+at),~{}~{}\psi_{3}(t)=\frac{at}{1+at},~{}~{}\psi_{4}(t)=\frac{-1}{1+at},$

where $a$ is positive number. In particular,

(1)

$\psi_{1}$ is a soft-thresholding penalty function HHM08 ; T96 . We have
$\psi^{\prime}_{1}(t)=1$ and $\psi^{\prime\prime}_{1}(t)=0$. Hence, it
satisfies Assumption (A).

(2)

$\psi_{2}$ is a logistic penalty function NNZC08 . We have

$\psi^{\prime}_{2}(t)=\frac{a}{1+at},~{}~{}\psi^{\prime\prime}_{2}(t)=-\frac{a^{2}}{(1+at)^{2}},$

which implies that $0<\lim_{t\to 0}\psi_{2}^{\prime}(t)=a$ and
$\lvert\psi^{\prime\prime}_{2}(t)\rvert\leq a^{2}$. Hence, it satisfies
Assumption (A).

(3)

$\psi_{3}$ is fraction penalty function CZh10 ; NNZC08 . We have

$\psi^{\prime}_{3}(t)=\frac{a}{(1+at)^{2}},~{}~{}\psi^{\prime\prime}_{3}(t)=-\frac{2a^{2}}{(1+at)^{3}},$

which implies that $0<\lim_{t\to 0}\psi_{3}^{\prime}(t)=a$ and
$\lvert\psi^{\prime\prime}_{3}(t)\rvert\leq 2a^{2}$. Hence, it satisfies
Assumption (A).

(4)

For function $\psi_{4}$, we have

$\psi^{\prime}_{4}(t)=\frac{a}{(1+at)^{2}},~{}~{}\psi^{\prime\prime}_{4}(t)=-\frac{2a^{2}}{(1+at)^{3}},$

which implies that $0<\lim_{t\to 0}\psi_{4}^{\prime}(t)=a$ and
$\lvert\psi^{\prime\prime}_{4}(t)\rvert\leq 2a^{2}$. Hence, it satisfies
Assumption (A).

##  Appendix B Proof of Lemma 3.1

In this appendix, we give a proof of Lemma 3.1.

Proof: By Assumption (B4), we know that there exists $\bar{d}\in\Re^{n+r}$
such that (26) and (27) hold. Meanwhile, we have from the formula of
$\Psi_{j}$ that

$\nabla_{(\omega,\lambda)}\Psi_{j}(\omega,\lambda)=\begin{cases}\nabla_{(\omega,\lambda)}\Phi_{j}(\omega,\lambda)&\text{if}~{}j\notin
I(\omega^{*})\\\ \left[\begin{array}[]{c}e_{j}\\\
0_{r}\end{array}\right]&\text{if}~{}j\in I(\omega^{*}),\end{cases}$ (80)

where $e_{j}$ is the $j$th standard unit vector in $\Re^{n}$ and $0_{r}$
denotes the zero vector in $\Re^{r}$. It is then clear from (80) and (26) that
$\bar{d}_{j}=0$ for all $j\in I(\omega^{*})$. Consequently, letting
$d\in\Re^{n-\lvert I(\omega^{*})\rvert+r}$ be the vector
$d:=(\bar{d})_{j\notin I(\omega^{*})}$, it follows from (26) and (27) that
equations (28) and (29) hold.

It remains to show that
$\\{\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})\\}_{j\notin
I(\omega^{*})}$ is linearly independent. To this end, note first that we have
from Assumption (B4) the linear independence of
$\\{\nabla_{(\omega,\lambda)}\Psi_{j}(\omega^{*},\lambda^{*})\\}_{j=1}^{n}$,
that is, the matrix

$M:=\left[\left(\nabla_{(\omega,\lambda)}\Psi_{j}(\omega^{*},\lambda^{*})\right)_{(j\notin
I(\omega^{*}))},\left(\nabla_{(\omega,\lambda)}\Psi_{j}(\omega^{*},\lambda^{*})\right)_{(j\in
I(\omega^{*}))}\right]\in\Re^{(n+r)\times n}$

has full column rank. Using (80) and switching the rows of $M$ so that the
first $\lvert I(\omega^{*})\rvert$ rows correspond to the index set
$I(\omega^{*})$, we have that the matrix

$\begin{bmatrix}(\nabla_{\hat{\omega}}\Phi_{j}(\omega^{*},\lambda^{*}))_{j\in
I(\omega^{*})}&E_{\lvert I(\omega^{*})\rvert}\\\
(\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*}))_{j\notin
I(\omega^{*})}&O_{(n-\lvert I(\omega^{*})\rvert+r)\times\lvert
I(\omega^{*})\rvert}\\\ \end{bmatrix}$

has full column rank, where $\hat{\omega}:=(\omega_{j})_{j\in I(\omega^{*})}$,
$E_{s}$ denotes the identity matrix of order $s$, and $O_{s\times t}$ is the
zero matrix of size $s\times t$. Since the upper and lower right blocks of the
above matrix are the identity matrix and zero matrix, respectively, a series
of elementary column operations leads us to conclude that

$\begin{bmatrix}O_{\lvert I(\omega^{*})\rvert\times(n-\lvert
I(\omega^{*})\rvert)}&E_{\lvert I(\omega^{*})\rvert}\\\
(\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*}))_{j\notin
I(\omega^{*})}&O_{(n-\lvert I(\omega^{*})\rvert+r)\times\lvert
I(\omega^{*})\rvert}\\\ \end{bmatrix}$

also has full column rank. As a consequence,
$\\{\nabla_{(\bar{\omega},\lambda)}\Phi_{j}(\omega^{*},\lambda^{*})\\}_{j\notin
I(\omega^{*})}$ is linearly independent, as desired.

## References

  * (1) Bergstra, J., Bengio, Y.: Random search for hyper-parameter optimization. J. Mach. Learn. Res. 13(10), 281–305 (2012)
  * (2) Bennett, K. P., Hu, J., Ji, X., Kunapuli, G., Pang, J-S.: Model selection via bilevel optimization. In The 2006 IEEE International Joint Conference on Neural Network Proceedings. 1922–1929 (2006)
  * (3) Bennett, K. P., Hu, J., Ji, X., Kunapuli, G., Pang, J-S.: Bilevel optimization and machine learning. Computational Intelligence: Research Frontiers. WCCI 2008. Lecture Notes in Computer Science, vol. 5050, Berlin, Heidelberg, 2008. Springer.
  * (4) Bian, W., Chen, X.: Worst-case complexity of smoothing quadratic regularization methods for non-Lipschitzian optimization. SIAM J. Optim. 23(3), 1718–1741 (2013)
  * (5) Bracken, J., McGill, J.: Mathematical programs with optimization problems in the constraints. Oper. Res. 21(1), 37–44 (1973)
  * (6) Chen, X.: Smoothing methods for nonsmooth, nonconvex minimization. Math. Program. 134(1), 71–99 (2012)
  * (7) Chen, C., Mangasarian, O. L.: A class of smoothing functions for nonlinear and mixed complementarity problems. Comput. Optim. Appl. 5(2), 97–138 (1996)
  * (8) Chen, X., Niu, L., Yuan, Y.: Optimality conditions and a smoothing trust region Newton method for nonLipschitz optimization. SIAM J. Optim. 23(3), 1528–1552 (2013)
  * (9) Chen, X., Xu, F., Ye, Y.: Lower bound theory of nonzero entries in solutions of $\ell_{2}$-$\ell_{p}$ minimization. SIAM J. Sci. Comput. 32(5), 2832–2852 (2010)
  * (10) Chen, X., Zhou, W.: Smoothing nonlinear conjugate gradient method for image restoration using nonsmooth nonconvex minimization. SIAM J. Imaging Sciences. 3(4), 765–790 (2010)
  * (11) Clarke, F. H.: Generalized gradients and applications. Trans. Amer. Math. Soc. 205, 247–262 (1975)
  * (12) Clarke, F. H.: Optimization and nonsmooth analysis. Classics Appl. Math. 5, SIAM, Philadelphia (1990)
  * (13) Colson, B., Marcotte, P., Savard, G.: An overview of bilevel optimization. Ann. Oper. Res. 153, 234–256 (2007)
  * (14) Dempe, S.: Annotated bibliography on bilevel programming and mathematical programs with equilibrium constraints. Optimization. 52(3), 333–359 (2003)
  * (15) Kruger, A. Y.: On Fréchet subdifferentials. J. Math. Sci. 116(3), 3325–3358 (2003)
  * (16) Kunisch, K., Pock, T.: A bilevel optimization approach for parameter learning in variational models. SIAM J. Imaging Sciences. 6(2), 938–983 (2013)
  * (17) Huang, J., Horowitz, J. L., Ma, S.: Asymptotic properties of bridge estimators in sparse high-dimensional regression models. Ann. Statist. 36(2) 587–613 (2008)
  * (18) Moore, G. M., Bergeron, C., Bennett, K. P.: Nonsmooth bilevel programming for hyperparameter selection. In IEEE International Conference on Data Mining Workshops. pp. 374–381 (2009)
  * (19) Moore, G. M., Bergeron, C., Bennett, K. P.: Model selection for primal SVM. 85, 175–208 (2011)
  * (20) Nikolova, M., Ng, M. K., Zhang, S., Ching, W.-K.: Efficient reconstruction of piecewise constant images using nonsmooth nonconvex minimization. SIAM J. Imaging Sciences. 1(1), 2–25 (2008)
  * (21) Ochs, P., Ranftl, R., Brox, T., Pock, T.: Bilevel optimization with nonsmooth lower level problems. In International Conference on Scale Space and Variational Methods in Computer Vision (SSVM) (2015)
  * (22) Ochs, P., Ranftl, R., Brox, T., Pock, T.: Techniques for gradient-based bilevel optimization with non-smooth lower level problems. J. Math. Imaging Vis. 56(2), 175–194 (2016)
  * (23) Okuno, T., Takeda, A., Kawana, A., Watanabe, M.: On $\ell_{p}$-hyperparameter learning via bilevel nonsmooth optimization. arXiv preprint arXiv:1806.01520v3 (2021), to appear in Journal of Machine Learning Research
  * (24) Rockafellar, R. T., Wets, R. J.-B.: Variational analysis. vol. 317. Springer Science & Business Media (2009)
  * (25) Saheya, B., Nguyen, C. T., Chen, J.-S.: Neural network based on systematically generated smoothing functions for absolute value equation. J. Appl. Math. Comput. 61, 533-558 (2019)
  * (26) Sinha, A., Malo, P., Deb, K.: A review on bilevel optimization: from classical to evolutionary approaches and applications. IEEE Trans. Evol. Comput. 22(2), 276–295 (2018)
  * (27) Tibshirani, R.: Regression shrinkage and selection via the Lasso. J. R. Statist. Soc. Ser. B. 58(1), 267–288 (1996)
  * (28) Tso, W. W., Burnak, B., Pistikopoulos, E. N.: HY-POP: Hyperparameter optimization of machine learning models through parametric programming. Comput. Chem. Eng. 139, 106902 (2020)
  * (29) Voronin, S., Ozkaya, G., Yoshida, D.: Convolution based smooth approximations to the absolute value function with application to non-smooth regularization. arXiv preprint arXiv:1408.6795v2 (2015)
  * (30) Wang, Z., Hutter, F., Zoghi, M., Matheson, D., de Freitas, N.: Bayesian optimization in a billion dimensions via random embeddings. J. Artif. Intell. Res. 55, 361–387 (2016)