# A proximal-proximal majorization-minimization algorithm for nonconvex
tuning-free robust regression problems

Peipei Tang, Chengjing Wang and Bo Jiang School of Computer and Computing
Science, Zhejiang University City College, Hangzhou 310015, China (Email:
tangpp@zucc.edu.cn). This author’s research is supported by the Natural
Science Foundation of Zhejiang Province of China under Grant No. LY19A010028,
the Zhejiang Science and Technology Plan Project of China (No. 2020C03091, No.
2021C01164) and the Scientific Research Foundation of Zhejiang University City
College (No. X-202112).Corresponding author, School of Mathematics, Southwest
Jiaotong University, No.999, Xian Road, West Park, High-tech Zone, Chengdu
611756, China (Email: renascencewang@hotmail.com). School of Computer Science,
Zhejiang University, No. 38, Zheda Road, Hangzhou 310027, China (Email:
22021105@zju.edu.cn).

###### Abstract

In this paper, we introduce a proximal-proximal majorization-minimization
(PPMM) algorithm for nonconvex tuning-free robust regression problems. The
basic idea is to apply the proximal majorization-minimization algorithm to
solve the nonconvex problem with the inner subproblems solved by a sparse
semismooth Newton (SSN) method based proximal point algorithm (PPA). We must
emphasize that the main difficulty in the design of the algorithm lies in how
to overcome the singular difficulty of the inner subproblem. Furthermore, we
also prove that the PPMM algorithm converges to a d-stationary point. Due to
the Kurdyka-Łojasiewicz (KL) property of the problem, we present the
convergence rate of the PPMM algorithm. Numerical experiments demonstrate that
our proposed algorithm outperforms the existing state-of-the-art algorithms.

Keywords:  Tuning-free robust regression problems, nonconvex regularizer,
proximal-proximal majorization-minimization algorithm, proximal point
algorithm, semismooth Newton method

##  1 Introduction

During the last decade, the estimation of high dimensional sparse statistical
models has attracted tremendous interests, in which the number of variables is
probably larger than the sample size. Consider a high-dimensional linear
regression model

$\displaystyle b=X\dot{\beta}+\varepsilon,$

where $b\in\mathcal{R}^{n}$ is an $n$-dimensional response vector,
$X\in\mathcal{R}^{n\times p}$ is a design matrix,
$\varepsilon\in\mathcal{R}^{n}$ is a random error vector. Sparse learning
based methods, which aim to minimize the fitting errors along with some sparse
regularization terms, have received considerable attention in recent years due
to their good performance and interpretability. The $\ell_{1}$ regularized
least square regression that is called LASSO [39] requires solving the
following optimization problem

$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\Big{\\{}\frac{1}{2n}\|X\beta-b\|^{2}+\lambda\|\beta\|_{1}\Big{\\}},$

where $\|\cdot\|$ is the Euclidean norm in ${\cal R}^{n}$ and $\lambda$
denotes the tuning parameter which controls the complexity of the model. The
Lasso estimator can achieve a near-oracle performance if the error vector
$\varepsilon$ is normally distributed and some suitable design conditions
hold. The theory of the Lasso reveals that the tuning parameter $\lambda$
relies on the deviation of the noise in order to recover $\dot{\beta}$.
However, it is of great challenge to estimate the deviation for large-scale
problems, especially when the feature dimension $p$ is much larger than the
sample size $n$. The square-root Lasso estimator

$\displaystyle\min_{\beta\in{\cal
R}^{p}}\Big{\\{}\|X\beta-b\|+\lambda\|\beta\|_{1}\Big{\\}}$

was proposed in [5] to avoid the aforementioned challenge. It is known (see
e.g., [4, 10]) that the square-root Lasso estimator can achieve the minimax
optimal rate of convergence under some suitable conditions which are
independent of the unknown noise level. Another obstacle is that the Gaussian
or sub-Gaussian error assumption can hardly be satisfied for high-dimensional
microarray data, climate data, insurance claim data, e-commerce data and many
other applications due to the heavy-tailed errors, which can affect the choice
of $\lambda$ and result in misleading results if we apply the standard
procedures directly. Therefore, Wang et al. [40] studied the following
$\ell_{1}$ regularized tuning-free robust regression model

$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\left\\{\frac{1}{n(n-1)}\sum_{1\leq
i<j\leq
n}\left|(b_{i}-X_{i}\beta)-(b_{j}-X_{j}\beta)\right|+\lambda\|\beta\|_{1}\right\\},$
(1)

where $X_{i}$ is the $i$-th row of the matrix $X$ and $b_{i}$ is the $i$-th
component of the vector $b$. The loss in (1) was originally from the classical
nonparametric statistics (see e.g., [18]) and is equivalent to Jaeckel’s
dispersion function (see [21]) with Wilcoxon scores. For simplicity, we call
the problem (1) the tuning-free robust Lasso problem. It has been shown in
[40] that the model (1) is very close to the Lasso for normal random errors
and is robust with substantial efficiency gain for heavy-tailed errors.

As we know, it is computationally challenge to solve the problems of Lasso,
square-root Lasso and tuning-free robust Lasso, though they are all convex.
Many numerical algorithms have been proposed to solve the Lasso problem,
including the accelerated proximal gradient (APG) method [3], the interior-
point method (IPM) [22], the least angle regression (LARS) [12] and the
semismooth Newton augmented Lagrangian method (SSNAL) [24], and so on. Due to
the nonsmooth loss in the objective function, the authors in [27] applied the
alternating direction method of multipliers (ADMM) to solve the square-root
Lasso problem. In a recent work [38], an efficient proximal majorization-
minimization algorithm (PMM) has been proposed to solve the square-root
regression problems with nonconvex regularizers. As for the tuning-free robust
Lasso problem, the authors in [40] reformulated this problem to an
$(n(n-1)+2p)$-dimensional linear programming problem with $(n(n-1)+2p)$
constraints. The well known IPM, which is implemented in some commercial
solvers such as Gurobi and Mosek, is a powerful algorithm for solving large-
scale problems. When the feature dimension $p$ and/or the sample size $n$ is
very large, the computational cost of linear programming is unacceptable.

Although the $\ell_{1}$ regularizer has many attractive properties, for large
regression coefficients the shrinkage introduced by the Lasso tends to be
significantly biased toward 0 and leads to biased estimate. Many alternative
regularizers such as the smoothly clipped absolute deviation (SCAD)
regularizer [14, 15] and the minimax concave penalty (MCP) [41] were designed
to diminish this bias. All these regularization regression models have the so-
called oracle property for individual variables, which means that the
corresponding regularized estimator is equal to the least squares estimator
assuming the model is known with high probability under appropriate
conditions. By making a local linear approximation to the regularizers, Zou
and Li [42] applied the LARS algorithm to solve these optimization problems.
It is known that many of these nonconvex regularizers, which are surrogate
sparsity functions, can be expressed as the difference of two convex functions
(DC) [1, 23]. By the DC property of the objective function, we can naturally
apply a majorization-minimization algorithm to solve the nonconvex problem.

In this paper, our main contribution is in proposing an efficient algorithm
for the tuning-free robust regression problems with nonconvex regularizations,
which takes the following form

$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\left\\{g(\beta):=\frac{1}{n(n-1)}\sum_{1\leq
i<j\leq n}\left|(b_{i}-X_{i}\beta)-(b_{j}-X_{j}\beta)\right|+\lambda
q_{1}(\beta)-q_{2}(\beta)\right\\},$ (2)

where $q_{1}:\mathcal{R}^{p}\rightarrow\mathcal{R}$ is a convex function with
a strongly semismooth proximal mapping and
$q_{2}:\mathcal{R}^{p}\rightarrow\mathcal{R}$ is a convex smooth function (the
dependence of $q_{2}$ on $\lambda$ has been dropped here). Since the proximal
mappings of many commonly used functions such as the $\ell_{1}$ and $\ell_{2}$
functions are strongly semismooth [30], the assumption on $q_{1}$ is mild. For
the DC structure of the regularization function, we design a two stage
proximal-proximal majorization-minimization (PPMM) algorithm to solve the
optimization problem (2). For each stage of the PPMM algorithm, the
corresponding subproblem is convex and contains two nonsmooth terms, i.e., the
loss function and the regularization function $q_{1}$. We apply a semismooth
Newton (SSN) method based proximal point algorithm (PPA) to solve the dual
problem of each subproblem. Note that the Clarke generalized Jacobian of the
subproblem is singular, and using the SSN method directly is not feasible. To
overcome this difficulty, we have to apply the SSN method to a sequentially
regularized subproblems.

The remaining parts of this paper are organized as follows. In Section 2, we
introduce some basic knowledge which will be used in this paper. In Section 3,
we briefly introduce the preconditioned PPA. In Section 4, we present the
details of the PPMM algorithm, and we also prove the convergence and the
convergence rate of the proposed algorithm based on the Kurdyka-Łojasiewicz
(KL) property. Furthermore, we prove that the algorithm converges to a
d-stationary point. In Section 5, we compare our algorithm with the existing
algorithms to demonstrate that the proposed algorithm can solve the tuning-
free robust regression problem efficiently.

###  1.1 Additional notations

Let $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ be the standard Euclidean
inner product and norm in the space $\mathcal{R}^{n}$. For any
$x\in\mathcal{R}^{n}$ and a given self-adjoint positive semidefinite matrix
$M\in\mathcal{R}^{n\times n}$, define $\|x\|_{M}:=\sqrt{\langle x,Mx\rangle}$.
The largest eigenvalue of $M$ is denoted by $\lambda_{\max}(M)$. Given a set
$C\subseteq\mathcal{R}^{n}$, the weighted distance of $x$ to $C$ is defined by
$\mbox{dist}_{M}(x,C):=\inf_{y\in C}\\{\|y-x\|_{M}\\}$. If $C=\emptyset$, we
have that $\mbox{dist}_{M}(x,C)=+\infty$ for all $x\in\mathcal{R}^{n}$. Denote
the identity matrix of order $n$ by $I_{n}$. If $M$ is an identity matrix, we
just omit the subscript matrix $M$. We define the indicator function
$\delta_{C}$ of the set $C$ by $\delta_{C}(x)=0$ if $x\in C$, otherwise
$\delta_{C}(x)=+\infty$. The intersection of all the convex sets containing
$C$ is called the convex hull of $C$ and is denoted by conv$(C)$.

##  2 Preliminaries

In this section, we present some basic preliminaries that will be used in this
paper.

For any real valued function
$r:\mathcal{R}^{n}\rightarrow\overline{\mathcal{R}}$, the conjugate to $r$ is
defined by

$\displaystyle r^{*}(x):=\sup_{v\in\textrm{dom}(r)}\Big{\\{}\langle
x,v\rangle-r(v)\Big{\\}}.$

We call $r$ a proper function if there exists at least one
$x\in\mathcal{R}^{n}$ such that $r(x)<+\infty$ and $r(x)>-\infty$ for all
$x\in\mathcal{R}^{n}$, or in other words, if $\mbox{dom}(r)$ is a nonempty set
and $r$ is finite; otherwise, it is improper. For a proper lower
semicontinuous function $r$, the Moreau envelope function $e_{\sigma r}$ and
the proximal mapping $\mbox{Prox}_{\sigma r}$ with parameter $\sigma>0$ of $r$
are defined respectively by

$\displaystyle e_{\sigma r}(x)$ $\displaystyle:=$
$\displaystyle\inf_{w\in\mathcal{R}^{n}}\Big{\\{}r(w)+\frac{1}{2\sigma}\|w-x\|^{2}\Big{\\}},$
$\displaystyle\mbox{Prox}_{\sigma r}(x)$ $\displaystyle:=$
$\displaystyle\mathop{\mbox{argmin}}_{w\in\mathcal{R}^{n}}\Big{\\{}r(w)+\frac{1}{2\sigma}\|w-x\|^{2}\Big{\\}}.$

When $r$ is also convex, it is known from Theorem 2.26 of [37] that the Moreau
envelope function $e_{\sigma r}$ is convex and continuously differentiable
with

$\displaystyle\nabla e_{\sigma
r}(x)=\frac{1}{\sigma}\big{(}x-\mbox{Prox}_{\sigma r}(x)\big{)}$

and the proximal mapping $\mbox{Prox}_{\sigma r}$ is single-valued and
continuous with the following Moreau’s identity (see, e.g., Theorem 31.5 of
[34]) holds

$\displaystyle\mbox{Prox}_{\sigma
r}(x)+\sigma\mbox{Prox}_{r^{*}/\sigma}(x/\sigma)=x,\quad\forall\
x\in\mathcal{R}^{n}.$

Consider a proper function $r:\mathcal{R}^{n}\rightarrow(-\infty,+\infty]$
with $x\in$ $\mbox{dom}(r)$, the regular subdifferential of $r$ at $x$ is
defined by

$\displaystyle\widehat{\partial}r(x):=\Big{\\{}v\in\mathcal{R}^{n}\ \Big{|}\
r(x^{\prime})\geq r(x)+\langle
v,x^{\prime}-x\rangle+o(\|x^{\prime}-x\|)\Big{\\}},$

which can be written equivalently as

$\displaystyle\widehat{\partial}r(x):=\Big{\\{}v\in\mathcal{R}^{n}\ \Big{|}\
\mathop{\liminf}_{x\neq x^{\prime}\rightarrow
x}\frac{r(x^{\prime})-r(x)-\langle
v,x^{\prime}-x\rangle}{\|x^{\prime}-x\|}\geq 0\Big{\\}}.$

If $x\notin\mbox{dom}(r)$, $\widehat{\partial}r(x)=\emptyset$. The limiting
subdifferential of $r$ at $x$ is defined by

$\displaystyle\partial r(x):=\Big{\\{}v\in\mathcal{R}^{n}\ \Big{|}\ \exists\
x^{\upsilon}\rightarrow x\ \mbox{and}\
v^{\upsilon}\in\widehat{\partial}r(x^{\upsilon})\ \mbox{with}\
v^{\upsilon}\rightarrow v\Big{\\}}.$

It is known from Theorem 8.6 of [37] that for $x\in\mbox{dom}(r)$,
$\widehat{\partial}r(x)$ and $\partial r(x)$ are closed with
$\widehat{\partial}r(x)$ convex and $\widehat{\partial}r(x)\subseteq\partial
r(x)$. If $r$ is a proper convex function with $x\in\mbox{dom}(r)$, the
regular subdifferential and the limiting subdifferential coincide with the set
of subgradients of $r$ at $x$ in the sense of convex analysis.

A function $R:\mathcal{R}^{n}\rightarrow\mathcal{R}^{m}$ is called
directionally differentiable at a point $x\in\mathcal{R}^{n}$ in a direction
$v\in\mathcal{R}^{n}$ if the limit

$\displaystyle R^{\prime}(x,v):=\lim_{t\downarrow 0}\frac{R(x+tv)-R(x)}{t}$

exists. If $R$ is directionally differentiable at $x$ in every direction
$v\in\mathcal{R}^{n}$, we say that $R$ is directionally differentiable at $x$.
If $m=1$, $R^{\prime}(x;v)$ coincides with the classic directional derivative
of a real valued function at $x$ in a direction $v$.

A point $x\in\mathcal{R}^{n}$ is called a directional-stationary
(d-stationary) point of a proper function $r$ if $r^{\prime}(x;v)\geq 0$ holds
for any $v\in\mathcal{R}^{n}$. According to generalized Fermat’s rule (see
e.g., Theorem 10.1 of [37]), the condition $0\in\widehat{\partial}r(x)$ is a
necessary condition for $x$ to be a local minimizer of $r$. Furthermore, if
$r$ is directionally differentiable at $x$ and locally Lipschitz continuous
near $x$, $0\in\widehat{\partial}r(x)$ if and only if $x$ is a d-stationary
point of $r$. If the function $r$ is proper and convex, the condition
$0\in\widehat{\partial}r(x)$ is not just necessary for a local minimum but
also sufficient for a global minimum.

Let $\mathcal{U}\subseteq\mathcal{R}^{n}$ be open and
$R:\mathcal{U}\rightarrow\mathcal{R}^{m}$ be a given vector-valued function
that is locally Lipschitz continuous. Let $\mathcal{U}_{R}$ be the subset of
$\mathcal{U}$ consisting of the points where $R$ is F(réchet)-differentiable
and $JR(x)\in\mathcal{R}^{m\times n}$ be the Jacobian matrix of $R$ at
$x\in\mathcal{U}_{R}$. It is known from Theorem 9.60 of [37] that the set
$\mathcal{U}\setminus\mathcal{U}_{R}$ is negligible. For any
$x\in\mathcal{U}$, define the B-subdifferential of $R$ at $x$ by

$\displaystyle\partial_{B}R(x):=\Big{\\{}V\in\mathcal{R}^{m\times n}\ \Big{|}\
\exists\ x^{\upsilon}\rightarrow x\ \mbox{with}\
x^{\upsilon}\in\mathcal{U}_{R}\ \mbox{and}\ JG(x^{\upsilon})\rightarrow
V\Big{\\}}.$

The Clarke subdifferential of $R$ at $x$ is defined by $\partial
R(x):=\mbox{conv}(\partial_{B}R(x))$.

For further discussion, we introduce the concept of semismoothness below. See
[31, 32] for more details.

######  Definition 2.1.

Given a locally Lipschitz continuous function
$R:\mathcal{U}\subseteq\mathcal{R}^{n}\rightarrow\mathcal{R}^{m}$ with an open
set $\mathcal{U}$ and a nonempty, compact valued and upper-semicontinuous
multifunction $\mathcal{K}:\mathcal{U}\rightrightarrows\mathcal{R}^{m\times
n}$, we call $R$ a semismooth function at $x\in\mathcal{U}$ with respect to
the multifunction $\mathcal{K}$ if

  1. (i)

$R$ is directionally differentiable at $x$;

  2. (ii)

$R(x+\Delta x)-R(x)-\Gamma\Delta x=o(\|\Delta x\|)$, for any
$\Gamma\in\mathcal{K}(x+\Delta x)$, $\Delta x\in\mathcal{R}^{n}$ and $\Delta
x\rightarrow 0$.

Furthermore, if (ii) is replaced by

$\displaystyle R(x+\Delta x)-R(x)-\Gamma\Delta x=O(\|\Delta x\|^{1+\gamma}),\
\mbox{for any}\ \Gamma\in\mathcal{K}(x+\Delta x),\ \Delta x\in\mathcal{R}^{n}\
\mbox{and}\ \Delta x\rightarrow 0,$

where $\gamma$ is a constant, then we call $R$ a $\gamma$-order (strongly if
$\gamma=1$) semismooth function at $x$ with respect to $\mathcal{K}$. The
function $R$ is said to be a semismooth function on $\mathcal{U}$ with respect
to $\mathcal{K}$ if it is semismooth everywhere on $\mathcal{U}$ with respect
to $\mathcal{K}$.

The Kurdyka-Łojasiewicz (KL) property (see e.g., [7, 8]) plays a central role
in our further analysis.

######  Definition 2.2.

A proper lower semicontinuous function
$r:\mathcal{R}^{n}\rightarrow(-\infty,+\infty]$ is said to have the KL
property at $x\in\mbox{dom}(\partial r)$ if there exist $\eta\in(0,+\infty]$,
a neighbour $\mathcal{U}$ of $x$ and a continuous concave function
$\varphi:[0,\eta)\rightarrow[0,+\infty)$ satisfying

  1. (1)

$\varphi(0)=0$;

  2. (2)

$\varphi$ is continuous at 0 and continuously differentiable on $(0,\eta)$;

  3. (3)

$\varphi^{\prime}(s)>0$, for all $0<s<\eta$

such that the KL inequality
$\varphi^{\prime}(r(x^{\prime})-r(x))\mbox{dist}(0,\partial r(x^{\prime}))\geq
1$ holds for any $x^{\prime}\in\mathcal{U}$ and
$r(x)<r(x^{\prime})<r(x)+\eta$. If $r$ satisfies the KL property at each point
of $\mbox{dom}(\partial r)$, then $r$ is called a KL function.

##  3 The preconditioned PPA

The PPA is a classical optimization algorithm, which dates back to [29] and
was extensively studied in [35, 36]. The preconditioned PPA which is a
generalization of the PPA was first studied by [26]. Consider a lower
semicontinuous proper convex function
$f:\mathcal{R}^{n}\rightarrow(-\infty,+\infty]$. Given a sequence of positive
real numbers $\\{\sigma_{k}\\}$ such that
$0<\sigma_{k}\uparrow\sigma_{\infty}\leq+\infty$, $x^{0}\in\mathcal{R}^{n}$
and $\\{M_{k}\\}$ being a sequence of self-adjoint positive definite matrix in
$\mathcal{R}^{n\times n}$ satisfying

$\displaystyle(1+\varrho_{k})M_{k}\succeq M_{k+1},\quad
M_{k}\succeq\lambda_{\min}I_{n},\ \forall\ k\geq
0,\quad\limsup_{k\rightarrow\infty}\lambda_{\max}(M_{k})=\lambda_{\infty}$

with $\\{\varrho_{k}\\}$ a nonnegative summable sequence and
$0<\lambda_{\min}\leq\lambda_{\infty}<+\infty$, the preconditioned PPA
generates a sequence $\\{x^{k}\\}$ approximately by

$\displaystyle x^{k+1}\approx\bar{x}^{k+1}=\mathop{\rm
argmin}_{x\in\mathcal{R}^{n}}\left\\{f(x)+\frac{1}{2\sigma_{k}}\|x-x^{k}\|_{M_{k}}^{2}\right\\}.$

If $M_{k}\equiv I_{n}$ for all $k\geq 0$, the preconditioned PPA becomes to
the classical PPA. There are two general criteria for the approximate
computation of $x^{k+1}$,

$\displaystyle(A)\quad\|x^{k+1}-\bar{x}^{k+1}\|_{M_{k}}\leq\varepsilon_{k},\quad
0\leq\varepsilon_{k},\quad\sum_{k=0}^{\infty}\varepsilon_{k}<+\infty,$
$\displaystyle(B)\quad\|x^{k+1}-\bar{x}^{k+1}\|_{M_{k}}\leq\delta_{k}\|x^{k+1}-x^{k}\|_{M_{k}},\quad
0\leq\delta_{k}<1,\quad\sum_{k=0}^{\infty}\delta_{k}<+\infty.$

A multifunction $\mathcal{F}:\mathcal{R}^{d}\rightrightarrows\mathcal{R}^{d}$
is said to be locally upper Lipschitz continuous at $x$ with modulus $\kappa$
which is independent of $x$ if there exists a neighbourhood $\mathcal{U}$ of
$x$ such that
$\mathcal{F}(y)\subseteq\mathcal{F}(x)+\kappa\|y-x\|\mathbb{B}_{d}$, $\forall\
y\in\mathcal{U}$, where $\mathbb{B}_{d}$ is a unit ball in $\mathcal{R}^{d}$.
A multifunction $\mathcal{F}$ is said to be piecewise polyhedral if its graph
$\mbox{gph}\mathcal{F}:=\\{(x,y)\ |\ y\in\mathcal{F}(x)\\}$ is the union of
finitely many polyhedral convex sets. The inverse of a piecewise polyhedral
multifunction is also piecewise polyhedral. It has been shown by [33] that a
piecewise polyhedral set-valued mapping has a fundamental locally upper
Lipschitz continuous property.

######  Proposition 3.1.

Let $\mathcal{F}:\mathcal{R}^{d}\rightrightarrows\mathcal{R}^{d}$ be a
piecewise polyhedral set-valued mapping. For any $x\in\mathcal{R}^{d}$, there
exists a constant $\kappa$ independent of $x$ such that $\mathcal{F}$ is
locally upper Lipschitz continuous with modulus $\kappa$.

We review some convergence results of the preconditioned PPA; one may also
refer to [26].

######  Theorem 3.1.

Suppose that $\Omega:=\\{x\ |\ 0\in\partial f(x)\\}\neq\emptyset$. Let
$\\{x^{k}\\}$ be any sequence generated by the preconditioned PPA under the
criterion $(A)$. Then the sequence $\\{x^{k}\\}$ is bounded and

$\displaystyle\mbox{dist}_{M_{k+1}}(x^{k+1},\Omega)\leq(1+\varrho_{k})\mbox{dist}_{M_{k}}(x^{k},\Omega)+(1+\varrho_{k})\varepsilon_{k},\quad\forall\
k\geq 0.$

In addition, the sequence $\\{x^{k}\\}$ converges to a point $x^{\infty}$ with
$0\in\partial f(x^{\infty})$.

######  Assumption 3.1.

The operator $\partial f$ satisfies the following error bound condition. For
any $\delta>0$, there exists $\kappa>0$ such that

$\displaystyle\mbox{dist}(x,(\partial
f)^{-1}(0))\leq\kappa\mbox{dist}(0,\partial
f(x)),\quad\forall\,x\in\\{x\,|\,\mbox{dist}(x,(\partial
f)^{-1}(0))\leq\delta\\}.$

######  Remark 3.1.

The error bound condition is critical to obtain the convergence rate of the
preconditioned PPA. It has been proven by Lemma 2.4 of [26] that a
multifunction satisfies the error bound condition if its inverse is locally
upper Lipschitz continuous at the origin.

######  Theorem 3.2.

Assume that $\Omega\neq\emptyset$ and Assumption 3.1 holds. Let $t$ be a
positive number satisfying
$t>\sum\limits_{k=0}^{\infty}\varepsilon_{k}(1+\varrho_{k})$ and the initial
point $x^{0}$ be a point such that

$\displaystyle\mbox{dist}_{M_{0}}(x^{0},\Omega)\leq\frac{t-\sum\limits_{k=0}^{\infty}\varepsilon_{k}(1+\varrho_{k})}{\prod\limits_{k=0}^{\infty}(1+\varrho_{k})}.$

Let $\\{x^{k}\\}$ be any sequence generated by the preconditioned PPA under
the criteria $(A)$ and $(B)$ with $\\{\sigma_{k}\\}$ nondecreasing. Then we
have

$\displaystyle\mbox{dist}_{M_{k+1}}(x^{k+1},\Omega)\leq\theta_{k}\mbox{dist}_{M_{k}}(x^{k},\Omega),\
\mbox{for all}\ k\geq 0,$

where

$\displaystyle\theta_{k}=(1+\varrho_{k})(1-\delta_{k})^{-1}\left(\delta_{k}+\frac{(1+\delta_{k})\kappa\lambda_{max}(M_{k})}{\sqrt{\sigma_{k}^{2}+\kappa^{2}\lambda_{max}^{2}(M_{k})}}\right)$

and

$\displaystyle\limsup_{k\rightarrow\infty}\theta_{k}=\theta_{\infty}=\frac{\kappa\lambda_{\infty}}{\sqrt{\sigma_{\infty}^{2}+\kappa^{2}\lambda_{\infty}^{2}}}<1\quad\mbox{with}\
\theta_{\infty}=0\ \mbox{if}\ \sigma_{\infty}=\infty.$

##  4 The PPMM algorithm

In this section, we introduce the PPMM algorithm to solve the nonconvex
tuning-free robust regression problem. Firstly, we obtain an initial point to
warm start the PPMM algorithm by solving a convex relaxation problem with
omitting the term $-q_{2}$ of the problem (2). We solve this subproblem by the
SSN method based PPA. Secondly, the solution of the nonconvex problem is
obtained by solving a series of majorizied problems. Specifically, the
majorized problem is constructed by linearizing the concave term $-q_{2}$ at
the current iterate point $\tilde{\beta}$. It is essentially important to
solve those subproblems efficiently and accurately. However, it is not easy to
accomplish this task. The main difficulty is that the Clarke generalized
Jacobian of the subproblem may be probably singular. To overcome this
difficulty, we introduce a special proximal term
$\frac{\tau}{2}\|\beta-\tilde{\beta}\|^{2}+\frac{\iota}{2}\|X(\beta-\tilde{\beta})\|^{2}$
($\tau>0$, $\iota>0$) and sequentially solve these regularized problems. After
introducing the proximal term, the dual of the problem can be written
explicitly as an unconstrained convex problem that can be solved by the SSN
method.

For $\tilde{\beta},\tilde{w}\in\mathcal{R}^{p}$, $\tau>0$ and $\iota>0$,
denote

$\displaystyle f(\beta;\tau,\iota,\tilde{\beta},\tilde{w})$ $\displaystyle:=$
$\displaystyle\frac{1}{n(n-1)}\sum_{1\leq i<j\leq
n}\left|(b_{i}-X_{i}\beta)-(b_{j}-X_{j}\beta)\right|+\lambda q_{1}(\beta)$ (3)
$\displaystyle-
q_{2}(\tilde{\beta})-\langle\tilde{w},\beta-\tilde{\beta}\rangle+\frac{\tau}{2}\|\beta-\tilde{\beta}\|^{2}+\frac{\iota}{2}\|X(\beta-\tilde{\beta})\|^{2}.$

Now we formally present the PPMM algorithm for the tuning-free robust
regression problem.

Algorithm 1 (PPMM):  Stage 1: (The preconditioned PPA with $M_{k}\equiv
I_{p}+X^{T}X$) Given $\tau_{1,0}>0$, $\beta^{0,0}\in\mathcal{R}^{p}$, set
$k=0$, iterate: Step 1.1. Apply Algorithm 2 to find an approximate solution
$u^{0,k+1}$ of the dual problem of the following optimization problem
$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\Big{\\{}f(\beta;\tau_{1,k},\tau_{1,k},\beta^{0,k},0)\Big{\\}}.$
Then $\beta^{0,k+1}=\mbox{Prox}_{\tau_{1,k}^{-1}\lambda
q_{1}}(-\tau_{1,k}^{-1}X^{T}u^{0,k+1}+\beta^{0,k})$. If a desired stopping
criterion is satisfied, let $\beta^{0}=\beta^{0,k+1}$, go to Stage 2;
otherwise go to Step 1.2. Step 1.2. Update
$\tau_{1,k+1}=\rho_{1,k}\tau_{1,k}$, $\rho_{1,k}\in(0,1)$. Set $k:=k+1$ and
return to Step 1.1. Stage 2: Given $\tau_{2,0}>0$ and $\iota_{2,0}>0$, set
$k=0$, iterate: Step 2.1. Apply Algorithm 2 to find an approximate solution
$u^{k+1}$ of the dual problem of the following optimization problem
$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\Big{\\{}f(\beta;\tau_{2,k},\iota_{2,k},\beta^{k},\nabla
q_{2}(\beta^{k}))\Big{\\}}.$ (4) Then
$\beta^{k+1}=\mbox{Prox}_{\tau_{2,k}^{-1}\lambda
q_{1}}(-\tau_{2,k}^{-1}X^{T}u^{k+1}+\tau_{2,k}^{-1}\nabla
q_{2}(\beta^{k})+\beta^{k})$. If a desired stopping criterion is satisfied,
terminate; otherwise go to Step 2.2. Step 2.2. Update
$\tau_{2,k+1}=\rho_{2,k}\tau_{2,k}$,
$\iota_{2,k+1}=\rho^{\prime}_{2,k}\iota_{2,k}$,
$\rho_{2,k}\in(0,1),\rho^{\prime}_{2,k}\in(0,1)$. Set $k:=k+1$ and return to
Step 2.1.

###  4.1 The subproblem of the PPMM algorithm

Given $\tau>0$, $\iota>0$ and $\tilde{\beta}\in\mathcal{R}^{p}$, consider the
following optimization problem

$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\Big{\\{}f(\beta;\tau,\iota,\tilde{\beta},\tilde{w})\Big{\\}}.$
(5)

Denote $h(y):=\frac{1}{n(n-1)}\sum\limits_{1\leq i<j\leq
n}\left|y_{i}-y_{j}\right|$, the problem (5) is equivalent to the following
minimization problem

$\displaystyle\min_{{\beta\in\mathcal{R}^{p}}\atop{y\in\mathcal{R}^{n}}}\Big{\\{}h(y)+\lambda
q_{1}(\beta)-\langle\tilde{w},\beta-\tilde{\beta}\rangle+\frac{\tau}{2}\|\beta-\tilde{\beta}\|^{2}+\frac{\iota}{2}\|y+b-X\tilde{\beta}\|^{2}\
\Big{|}\ X\beta-y-b=0\Big{\\}}.$ (6)

We can also write the dual problem of (6) as

$\displaystyle\min_{u\in\mathcal{R}^{n}}\Big{\\{}\phi(u):=\frac{\iota}{2}\|\iota^{-1}u+X\tilde{\beta}-b\|^{2}-\frac{\iota}{2}\|\iota^{-1}u+X\tilde{\beta}-b-\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u+X\tilde{\beta}-b)\|^{2}$
$\displaystyle\quad\quad\quad-h(\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u+X\tilde{\beta}-b))+\frac{\tau}{2}\|-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta}\|^{2}$
(7)
$\displaystyle\quad\quad\quad-\frac{\tau}{2}\|-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta}-\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta})\|^{2}$
$\displaystyle\quad\quad\quad-\lambda q_{1}(\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta}))+\langle
u,b\rangle\Big{\\}}.$

Let $\bar{u}$ be a solution of the problem (4.1), then we have
$\bar{\beta}=\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta})$ and
$\bar{y}=\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u+X\tilde{\beta}-b)$, which is a
solution of the primal problem (6).

The objective function of the problem (4.1) is smooth and finding the optimal
solution of the problem (4.1) is equivalent to solving the following system of
equations

$\displaystyle\nabla\phi(u)=\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u+X\tilde{\beta}-b)-X\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta})+b=0.$ (8)

One may apply the SSN method to solve the aforementioned problem (8). However,
the Clarke generalized Jacobian of the problem (8) may be probably singular
and the SSN method cannot be directly applied. To overcome this difficulty, we
deal with the problem (4.1) in an alternative way by applying the PPA with the
inner subproblem solved by the SSN method.

Given $u^{0}\in\mathcal{R}^{n}$, let $\\{u^{i}\\}$ be the sequence generated
by the PPA. We adopt

$\displaystyle(B^{\prime\prime})\quad\|\nabla\phi(u^{i+1})+\sigma_{i}^{-1}(u^{i+1}-u^{i})\|\leq\frac{\delta_{i}}{\sigma_{i}}\|u^{i+1}-u^{i}\|,\
0\leq\delta_{i}<1,\ \sum_{i=0}^{\infty}\delta_{i}<\infty$

as the stopping criterion of each iteration. We apply the SSN method to
compute an approximate solution of each subproblem. Let
$\hat{\phi}(u;\sigma_{i},u^{i}):=\phi(u)+\frac{1}{2\sigma_{i}}\|u-u^{i}\|^{2}$.
Since the mappings $\mbox{Prox}_{\iota^{-1}h}(\cdot)$ and
$\mbox{Prox}_{\tau^{-1}\lambda q_{1}}(\cdot)$ are both Lipschitz continuous,
the following multifunction

$\displaystyle\hat{\partial}^{2}\hat{\phi}(u;\sigma_{i},u^{i})$
$\displaystyle:=$
$\displaystyle\iota^{-1}\partial\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u+X\tilde{\beta}-b)$
$\displaystyle+\tau^{-1}X\partial\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta})X^{T}+\sigma_{i}^{-1}I_{n}$

is well defined. It is known from [20] that

$\displaystyle\partial^{2}\hat{\phi}(u;\sigma_{i},u^{i})(d)=\hat{\partial}^{2}\hat{\phi}(u;\sigma_{i},u^{i})(d),\quad\forall\
d\in\mathcal{R}^{n},$

where $\partial^{2}\hat{\phi}(u;\sigma_{i},u^{i})$ is the Clarke generalized
Jacobian of $\hat{\phi}(u;\sigma_{i},u^{i})$ at $u$. Let us select
$U\in\partial\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u+X\tilde{\beta}-b)$ and
$V\in\partial\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u+\tau^{-1}\tilde{w}+\tilde{\beta})$, we obtain
$H:=\iota^{-1}U+\tau^{-1}XVX^{T}+\sigma_{i}^{-1}I_{n}\in\partial^{2}\hat{\phi}(u;\sigma_{i},u^{i})$.

Now we introduce the PPA and the SSN method as below.

Algorithm 2 (PPA):  Given $\sigma_{0}>0$, choose $u^{0}\in\mathcal{R}^{n}$.
Set $i=0$ and iterate: Step 1. Apply Algorithm 3 to solve the following
problem $\displaystyle u^{i+1}\approx\mathop{\rm
argmin}_{u\in\mathcal{R}^{n}}\hat{\phi}(u;\sigma_{i},u^{i}).$ (9) Step 2. If a
desired stopping criterion is satisfied, terminate; otherwise, update
$\sigma_{i+1}=\rho_{i}\sigma_{i}$ with $\rho_{i}\geq 1$ and return to Step 1.

Algorithm 3 (SSN):  Given $\iota,\tau,\sigma>0,$
$\tilde{\beta},\tilde{w}\in\mathcal{R}^{p}$, $\mu\in(0,\frac{1}{2}),\
\overline{\eta}\in(0,1),\varsigma\in(0,1],\textrm{ and }\delta\in(0,1)$,
choose $u^{0}\in\mathcal{R}^{n}$. Set $j=0$ and iterate: Step 1. Choose
$U^{j}\in\partial\mbox{Prox}_{\iota^{-1}h}(\iota^{-1}u^{j}+X\tilde{\beta}-b)$
and $V^{j}\in\partial\mbox{Prox}_{\tau^{-1}\lambda
q_{1}}(-\tau^{-1}X^{T}u^{j}+\tau^{-1}\tilde{w}+\tilde{\beta})$. Let
$H^{j}=\iota^{-1}U^{j}+\tau^{-1}XV^{j}X^{T}+\sigma^{-1}I_{n}$ and find the
solution $u^{j}$ of the following linear system $\displaystyle H^{j}\Delta
u=-\nabla\hat{\phi}(u^{j};\sigma,u^{0})$ such that $\displaystyle\|H^{j}\Delta
u^{j}+\nabla\hat{\phi}(u^{j};\sigma,u^{0})\|\leq\eta_{j}:=\min(\overline{\eta},\|\nabla\hat{\phi}(u^{j};\sigma,u^{0})\|^{1+\varsigma}).$
(10) Step 2. Set $\alpha_{j}=\delta^{m_{j}}$, where $m_{j}$ is the first
nonnegative integer $m$ for which
$\displaystyle\hat{\phi}(u^{j}+\delta^{m}\Delta
u^{j};\sigma,u^{0})\leq\hat{\phi}(u^{j};\sigma,u^{0})+\mu\delta^{m}\langle\nabla\hat{\phi}(u^{j};\sigma,u^{0}),(\Delta
u^{j})\rangle.$ Step 3. Set $u^{j+1}=u^{j}+\alpha_{j}\Delta u^{j}$.

Since $\mbox{Prox}_{\iota^{-1}h}(\cdot)$ is a piecewise affine function, it is
strongly semismooth due to Proposition 7.4.7 of [13]. We have assumed that
$\mbox{Prox}_{q_{1}}(\cdot)$ is strongly semismooth. Therefore the function
$\nabla\hat{\phi}(\cdot;\sigma,u^{0})$ is strongly semismooth. Hence we can
establish the result about the local convergence rate of the SSN method. We
just present the main result here without proof.

######  Theorem 4.1.

Suppose that $\sigma<+\infty$. The sequence $\\{u^{j}\\}$ generated by the SSN
method converges to the unique solution $\bar{u}$ of
$\nabla\hat{\phi}(v;\sigma,u^{0})=0$ with

$\displaystyle\|u^{j+1}-\bar{u}\|=\mathcal{O}(\|u^{j}-\bar{u}\|^{1+\varsigma}).$

###  4.2 The computation of $\mbox{Prox}_{\iota^{-1}h}(x)$ and its generalized
Jacobian

In the implementation of the SSN method, it is essentially important to
compute the structure of $\partial\mbox{Prox}_{\iota^{-1}h}(x)$. In this
subsection, we introduce some basic results about the proximal mapping of $h$
and its corresponding generalized Jacobian which have been considered by [28].
Let $\rho:=\frac{1}{\iota n(n-1)}$, the proximal mapping related to
$\iota^{-1}h$ can be computed by the following proposition.

######  Proposition 4.1.

Let $x\in\mathcal{R}^{n}$ be an arbitrarily given vector and
$P_{x}\in\mathcal{R}^{n\times n}$ be a corresponding permutation matrix such
that $x^{\downarrow}=P_{x}x$ and $x^{\downarrow}_{1}\geq
x^{\downarrow}_{2}\geq\cdots\geq x^{\downarrow}_{n}$. Then the proximal
mapping related to $\iota^{-1}h$ can be computed by

$\displaystyle\mbox{Prox}_{\iota^{-1}h}(x)=P_{x}^{T}\Pi_{D}(x^{\downarrow}-\rho
w),$

where $w\in\mathcal{R}^{n}$ with $w_{i}=n-2i+1$, $i=1,2,\ldots,n$,
$D=\\{x\in\mathcal{R}^{n}\ |\ Bx\geq 0\\}$ with $B$ being a matrix such that
$Bx=[x_{1}-x_{2};\ldots;x_{n-1}-x_{n}]\in\mathcal{R}^{n-1}$ and
$\Pi_{D}(\cdot)$ is the metric projection onto $D$ which can be computed by
the pool-adjacent-violators algorithm (see [6] for more details).

Now we introduce the generalized Jacobian of
$\mbox{Prox}_{\iota^{-1}h}(\cdot)$, denoted by $\mathcal{M}(\cdot)$. It is
known that $\mbox{Prox}_{\iota^{-1}h}(\cdot)$ is strongly semismooth with
respect to $\mathcal{M}(\cdot)$. The details of $\mathcal{M}(\cdot)$ can be
found in [28].

For a given vector $y\in\mathcal{R}^{n}$, the active index set is denoted by

$\displaystyle\mathcal{I}_{A}(y):=\Big{\\{}i\ \Big{|}\ B_{i}\Pi_{D}(y)=0,\
i=1,2,\ldots,n-1\Big{\\}},$

where $B_{i}$ is the $i$-th row of the matrix $B$. Based on the active index
set $\mathcal{I}_{A}(x^{\downarrow}-\rho w)$, we define a diagonal matrix
$\Sigma\in\mathcal{R}^{(n-1)\times(n-1)}$ with each diagonal element taking
the following value

$\displaystyle\Sigma_{ii}=\left\\{\begin{array}[]{cc}1,&i\in\mathcal{I}_{A}(x^{\downarrow}-\rho
w),\\\ 0,&\mbox{otherwise,}\end{array}\right.\quad i=1,2,\ldots,n-1.$

Now we are ready to get the following proposition which proposes a way to find
a computational element of $\mathcal{M}(\cdot)$.

######  Proposition 4.2.

For a given vector $x\in\mathcal{R}^{n}$, we have

$\displaystyle P_{x}^{T}\big{(}I_{n}-B^{T}(\Sigma
BB^{T}\Sigma)^{\dagger}B\big{)}P_{x}\in\mathcal{M}(x),$

where $(\cdot)^{\dagger}$ denotes the pseduoinverse.

###  4.3 The convergence analysis of the PPMM algorithm

For Stage 1 of the PPMM algorithm, it solves a convex and polyhedral problem,
hence the convergence results follow from Theorem 3.1 and 3.2.

Now we focus on the convergence analysis of Stage 2 of the PPMM algorithm on
the nonconvex problem. Let $r^{k+1}$ be the residual vector of the gradient
for the objective function of the dual problem related to the problem (4) at
$u^{k+1}$. By the stopping criterion ($B^{\prime\prime}$) of the PPA, we know
that $r^{k+1}$ converges to zero as $\sigma_{i}\rightarrow\infty$. It follows
that $r^{k+1}$ satisfies

$\displaystyle
2h(r^{k+1})+\iota_{2,k}\|r^{k+1}\|^{2}\leq\frac{\iota_{2,k}}{2}\|X(\beta^{k+1}-\beta^{k})\|^{2},$
(12)

if $\sigma_{i}$ is large enough. For simplicity, denote

$\displaystyle f_{k}(\beta)$ $\displaystyle:=$ $\displaystyle
h(X\beta-b+r^{k+1})+\lambda q_{1}(\beta)-q_{2}(\beta^{k})-\langle\nabla
q_{2}(\beta^{k}),\beta-\beta^{k}\rangle$
$\displaystyle+\frac{\tau_{2,k}}{2}\|\beta-\beta^{k}\|^{2}+\frac{\iota_{2,k}}{2}\|X(\beta-\beta^{k})+r^{k+1}\|^{2}.$

Then we can rewrite Stage 2 of the PPMM algorithm as

$\beta^{k+1}=\mathop{\mbox{argmin}}\limits_{\beta\in\mathcal{R}^{p}}\Big{\\{}f_{k}(\beta)\Big{\\}}.$

We first introduce the following descent property of the series
$\\{g(\beta^{k})\\}$.

######  Lemma 4.1.

Let $\\{\beta^{k}\\}$ be a sequence generated by the PPMM algorithm. We have
the following descent property

$\displaystyle g(\beta^{k})\geq
g(\beta^{k+1})+\frac{\tau_{2,k}}{2}\|\beta^{k+1}-\beta^{k}\|^{2}.$

###### Proof.

Due to the following property of the function $h$

$\displaystyle|h(x)-h(y)|\leq h(x-y),\ \mbox{for}\ \forall\
x,y\in\mathcal{R}^{n},$ (13)

we have

$\displaystyle f_{k}(\beta^{k})\leq
g(\beta^{k})+h(r^{k+1})+\frac{\iota_{2,k}}{2}\|r^{k+1}\|^{2}.$

On the other hand, since $q_{2}$ is a convex function and $\beta^{k+1}$ is a
minimizer of the function $f_{k}$, we obtain

$\displaystyle f_{k}(\beta^{k})$ $\displaystyle\geq$ $\displaystyle
f_{k}(\beta^{k+1})$ $\displaystyle\geq$ $\displaystyle
h(X\beta^{k+1}-b+r^{k+1})+\lambda
q_{1}(\beta^{k+1})-q_{2}(\beta^{k+1})+\frac{\tau_{2,k}}{2}\|\beta^{k+1}-\beta^{k}\|^{2}$
$\displaystyle+\frac{\iota_{2,k}}{2}\|X(\beta^{k+1}-\beta^{k})+r^{k+1}\|^{2}$
$\displaystyle\geq$ $\displaystyle
g(\beta^{k+1})+\frac{\tau_{2,k}}{2}\|\beta^{k+1}-\beta^{k}\|^{2}+\frac{\iota_{2,k}}{2}\|X(\beta^{k+1}-\beta^{k})\|^{2}-h(r^{k+1})-\frac{\iota_{2,k}}{2}\|r^{k+1}\|^{2}.$

Combining those above two results, we derive

$\displaystyle g(\beta^{k})\geq
g(\beta^{k+1})+\frac{\tau_{2,k}}{2}\|\beta^{k+1}-\beta^{k}\|^{2}+\frac{\iota_{2,k}}{2}\|X(\beta^{k+1}-\beta^{k})\|^{2}-2h(r^{k+1})-\iota_{2,k}\|r^{k+1}\|^{2}.$

By the condition (12), the desired result follows. ∎

The following lemma is similar to Lemmas 5 and 6 in [9], and we present it
here without proof.

######  Lemma 4.2.

Let $\bar{\beta}\in\mathcal{R}^{p}$. Then $\bar{\beta}$ is a d-stationary
point of (2) if and only if there exist $\tau,\iota\geq 0$ such that
$\bar{\beta}$ is a solution of the following minimization problem

$\displaystyle\min_{\beta\in\mathcal{R}^{p}}\Big{\\{}f(\beta;\tau,\iota,\bar{\beta},\nabla
q_{2}(\bar{\beta}))\Big{\\}}.$

Now we introduce the convergence result of the PPMM algorithm.

######  Theorem 4.2.

Assume that the objective function in (2) is bounded below and
$\\{\tau_{2,k}\\}$ and $\\{\iota_{2,k}\\}$ are positive convergent sequences.
Then every cluster point of the sequence $\\{\beta^{k}\\}$ generated by the
PPMM algorithm is a d-stationary point of (2).

###### Proof.

It is known from Lemma 4.1 that $\\{g(\beta^{k})\\}$ is a non-increasing
sequence. Therefore, the sequences $\\{g(\beta^{k})\\}$ and
$\\{\|\beta^{k+1}-\beta^{k}\|\\}$ both converge and $\lim\limits_{k\rightarrow
0}\|\beta^{k+1}-\beta^{k}\|=0$ due to the assumption of the objective
function. Let $\\{\beta^{k}\\}_{k\in\mathcal{K}}$ be a subsequence of
$\\{\beta^{k}\\}$ with
$\lim\limits_{k(\in\mathcal{K})\rightarrow\infty}\beta^{k}=\beta^{\infty}$. We
will prove that $\beta^{\infty}$ is a d-stationary point of (2).

By the definition of $\beta^{k+1}$, we have

$\displaystyle f(\beta;\tau_{2,k},\iota_{2,k},\beta^{k},\nabla
q_{2}(\beta^{k}))+h(r^{k+1})+\frac{\iota_{2,k}}{2}\|r^{k+1}\|^{2}$
$\displaystyle\geq$ $\displaystyle f_{k}(\beta)\geq f_{k}(\beta^{k+1})$
$\displaystyle\geq$ $\displaystyle
f(\beta^{k+1};\tau_{2,k},\iota_{2,k},\beta^{k},\nabla
q_{2}(\beta^{k}))-h(r^{k+1})-\frac{\iota_{2,k}}{2}\|r^{k+1}\|^{2},\quad\forall\
\beta\in\mathcal{R}^{p}.$

Let $\tau_{\infty}=\lim\limits_{k\rightarrow\infty}\tau_{2,k}$ and
$\iota_{\infty}=\lim\limits_{k\rightarrow\infty}\iota_{2,k}$. Since
$\|r^{k}\|\rightarrow 0$ due to (12) and $h$ is continuous, by letting
$k(\in\mathcal{K})\rightarrow\infty$ we derive

$\displaystyle f(\beta;\tau_{\infty},\iota_{\infty},\beta^{\infty},\nabla
q_{2}(\beta^{\infty}))\geq
f(\beta^{\infty};\tau_{\infty},\iota_{\infty},\beta^{\infty},\nabla
q_{2}(\beta^{\infty})).$

Therefore,
$\beta^{\infty}\in\mathop{\mbox{argmin}}\limits_{\beta\in\mathcal{R}^{p}}\Big{\\{}f(\beta;\tau_{\infty},\iota_{\infty},\beta^{\infty},\nabla
q_{2}(\beta^{\infty}))\Big{\\}}$ and the desired result follows from Lemma
4.2. ∎

Based on an isolation assumption of the cluster point or the KL property
assumption, we establish the local convergence rate of the sequence
$\\{\beta^{k}\\}$ in the following theorem.

######  Theorem 4.3.

Assume that the objective function in (2) is bounded below. Let
$\\{\beta^{k}\\}$ be the sequence generated by the PPMM algorithm with
$\mathcal{C}^{\infty}$ as the set of all its cluster points. The whole
sequence $\\{\beta^{k}\\}$ will converge to an element of
$\mathcal{C}^{\infty}$, if one of the following conditions holds:

  1. (I)

The set $\mathcal{C}^{\infty}$ contains an isolated element;

  2. (II)

$\\{\beta^{k}\\}$ is a bounded sequence. The objective function $g$ in (2) has
the KL property at $\bar{\beta}$ and $\nabla q_{2}$ is locally Lipschitz
continuous near $\bar{\beta}$ for any $\bar{\beta}\in\mathcal{C}^{\infty}$.

Furthermore, based on the condition (2), let
$\lim\limits_{k\rightarrow\infty}\beta^{k}=\beta^{\infty}\in\mathcal{C}^{\infty}$
and the function $g$ satisfies the KL property at $\beta^{\infty}$ with an
exponent $\gamma\in[0,1)$, we have

  1. (a)

The sequence $\\{\beta^{k}\\}$ converges in a finite number of steps, if
$\gamma=0$;

  2. (b)

The sequence $\\{\beta^{k}\\}$ converges R-linearly, i.e., there exist $\mu>0$
and $\theta\in[0,1)$ such that
$\|\beta^{k}-\beta^{\infty}\|\leq\mu\theta^{{}^{k}}$, if
$0<\gamma\leq\frac{1}{2}$ and $k$ is sufficiently large;

  3. (c)

The sequence $\\{\beta^{k}\\}$ converges R-sublinearly, i.e., there exists
$\mu>0$ such that $\|\beta^{k}-\beta^{\infty}\|\leq\mu
k^{-\frac{1-\gamma}{2\gamma-1}}$, if $\frac{1}{2}<\gamma<1$ and $k$ is
sufficiently large.

###### Proof.

It is already known from the proof of Theorem 4.2 that
$\lim\limits_{k\rightarrow\infty}\|\beta^{k+1}-\beta^{k}\|=0$. Under the
condition (I), the sequence converges to the isolated element of
$\mathcal{C}^{\infty}$ due to Proposition 8.3.10 of [13]. In order to prove
the convergence rate of the sequence $\\{\beta^{k}\\}$ based on the condition
(II), we first describe some properties as below.

  1. (i)

The sequence $\\{g(\beta^{k})\\}$ is descent, i.e., $g(\beta^{k+1})\geq
g(\beta^{k})+\frac{\tau_{2,k}}{2}\|\beta^{k+1}-\beta^{k}\|^{2}$;

  2. (ii)

There exists a subsequence $\\{\beta^{k}\\}_{k\in\mathcal{K}}$ of
$\\{\beta^{k}\\}$ such that
$\lim\limits_{k(\in\mathcal{K})\rightarrow\infty}\beta^{k}=\beta^{\infty}$ and
$\lim\limits_{k(\in\mathcal{K})\rightarrow\infty}g(\beta^{k})=g(\beta^{\infty})$;

  3. (iii)

There exist a constant $K>0$ and $\varepsilon^{k+1}\in\partial g(\beta^{k+1})$
such that $\|\varepsilon^{k+1}\|\leq K\|\beta^{k+1}-\beta^{k}\|$, for $k$
sufficiently large.

The first two properties are already known. To establish the third property,
let $\bar{\varepsilon}^{k+1}=\nabla q_{2}(\beta^{k})-\nabla
q_{2}(\beta^{k+1})-\tau_{2,k}(\beta^{k+1}-\beta^{k})-\iota_{2,k}X^{T}(X(\beta^{k+1}-\beta^{k})+r^{k+1})$.
Then $\bar{\varepsilon}^{k+1}\in\partial g_{k}(\beta^{k+1})$ with
$g_{k}(\beta)=h(X\beta-b+r^{k+1})+\lambda q_{1}(\beta)-q_{2}(\beta)$. Since
the function $h$ is polyhedral, it is known from Theorem 23.8 and Theorem 23.9
of [34] that $\partial g_{k}(\beta^{k+1})=X^{T}\partial
h(X\beta^{k+1}-b+r^{k+1})+\partial q(\beta^{k+1})$, where $q(\beta):=\lambda
q_{1}(\beta)-q_{2}(\beta)$. Furthermore, $\partial h$ is a polyhedral
multifunction and locally upper Lipschitz continuous at the point
$X\beta^{k+1}-b$. Therefore, there exists $\kappa>0$ such that $\partial
h(X\beta^{k+1}-b+r^{k+1})\subseteq\partial
h(X\beta^{k+1}-b)+\kappa\|r^{k+1}\|\mathbb{B}_{p}$ with $\mathbb{B}_{p}$
denoting the unit ball in $\mathcal{R}^{p}$. We can find
$d\bar{\varepsilon}^{k+1}\in\mathcal{R}^{p}$ with
$\|d\bar{\varepsilon}^{k+1}\|\leq\kappa\|X\|\|r^{k+1}\|$ such that
$\varepsilon^{k+1}:=\bar{\varepsilon}^{k+1}+d\bar{\varepsilon}^{k+1}\in\partial
g(\beta^{k+1})$. Since $\nabla q_{2}$ is locally Lipschitz continuous at all
$\beta\in\mathcal{C}^{\infty}$ and
$\|r^{k+1}\|\leq\frac{1}{\sqrt{2}}\|X\|\|\beta^{k+1}-\beta^{k}\|$ according to
(12), we can find some $K>0$ such that $\|\varepsilon^{k+1}\|\leq
K\|\beta^{k+1}-\beta^{k}\|$. Combining properties (i)-(iii) and Proposition 4
of [7], the desired result follows easily.

∎

##  5 Numerical experiments

In this section, we implement some numerical experiments to demonstrate the
efficiency of our PPMM algorithm for the tuning-free robust regression
problems. The numerical experiments are implemented on a PC (Intel Core 2 Duo
2.6 GHz with 4 GB RAM) on two types of data sets. The first type of data set
is generated randomly in the low sample high dimension setting which is
similar to that in [40]. The synthetic data is from the regression model

$\displaystyle Y_{i}=X_{i}^{T}\dot{\beta}+\varsigma_{i},\ i=1,2,\ldots,n,$

where $X_{i}$ is generated from a $p$-dimensional multivariate normal
distribution $N_{p}(0,\Sigma)$ and is independent of $\varsigma_{i}$. In the
first six examples, we set
$\dot{\beta}=(\sqrt{3},\sqrt{3},\sqrt{3},\textbf{0}_{p-3})$ with
$\textbf{0}_{p-3}$ a $(p-3)$-dimensional vector of zeros. The covariance
matrix $\Sigma$ is symmetric with $\Sigma_{ij}=0.5$ for $i\neq j$ and
$\Sigma_{ij}=1$ for $i=j$. The noise $\varsigma_{i}$ is from the following
different distributions.

  1. (1)

The normal distribution with mean 0 and variance 0.25;

  2. (2)

The normal distribution with mean 0 and variance 1;

  3. (3)

The normal distribution with mean 0 and variance 2;

  4. (4)

The mixture normal distribution $0.95N(0,1)+0.05N(0,100)$;

  5. (5)

The $t$ distribution with 4 degree of freedom;

  6. (6)

The standard Cauchy distribution.

The settings in the next six examples are the same as that in the first six
examples except that
$\dot{\beta}=(2,2,2,2,1.75,1.75,1.75,1.5,1.5,1.5,1.25,1.25,1.25,1,1,1,0.75,0.75,0.75,0.5,$
$0.5,0.5,0.25,0.25,0.25,\textbf{0}_{p-25})$. We set the random seed of the
$i$-th example as $i$ ($i=1,2,\ldots,12$).

We also implement our experiments on some large-scale data sets obtained from
the KEEL-dataset repository [2]. As that in [24], we apply the method in [19]
and use the polynomial basis functions to expand the features of those data
sets with the last digit in the names of the data sets , e.g., baseball5,
concrete7, dee10 and puma32h51tst3 denoting the order of the polynomial used
to expand the features. Let $\hat{\beta}$ be the vector obtained by sorting
$\beta$ such that
$|\hat{\beta}_{1}|\geq|\hat{\beta}_{2}|\geq\cdots\geq|\hat{\beta}_{p}|$. In
our numerical experiments, we define the number of nonzero elements as the
minimal $k$ satisfying

$\displaystyle\sum_{i=1}^{k}|\hat{\beta}_{i}|\geq 0.9999\|\beta\|_{1}.$

For simplicity, we use “$s\ \mbox{sign}(t)|t|$” to denote a number of the form
“$s\times 10^{t}$”, e.g., 1.0-3 denotes $1.0\times 10^{-3}$.

###  5.1 Numerical experiments for the $\ell_{1}$ regularized tuning-free
robust regression problem

In the subsection, we implement some performances of the ADMM, Gurobi and our
proposed PPMM algorithm of Stage 1 for the $\ell_{1}$ regularized tuning-free
robust regression problem (1). We adopt the relative KKT residual

$\displaystyle\eta_{kkt}:=\max\left\\{\frac{\|y-\mbox{Prox}_{h}(u+y)\|}{1+\|y\|},\frac{\|x-\mbox{Prox}_{\lambda\|\cdot\|_{1}}(x-A^{T}u)\|}{1+\|x\|},\frac{\|Ax-
y-b\|}{1+\|y\|}\right\\}$

to measure the accuracy of our PPMM algorithm and the ADMM. Our PPMM algorithm
and the ADMM are terminated if the desired relative KKT residual
$\eta_{kkt}<\mbox{Tol}:=10^{-6}$ or the number of the iterations reaches the
maximum of $N_{\max}$. In our performances, we set $N_{max}$ as 200 for the
PPMM algorithm and 40000 for the ADMM. As for Gurobi, we use the barrier
algorithm (without presolve and crossover) and set the relevant tolerance
BarConvTol$=10^{-6}$. All of these performances are stopped if the running
time reaches the pre-set maximum of 4 hours. The tuning parameter $\lambda$ is
obtained by simulation based on 1000 repetitions with the same parameter
setting as that in [40].

####  5.1.1 ADMM

In this subsection, we introduce the implementation details of the ADMM for
the problem (1). We first introduce auxiliary variables and reformulate (1) as
follows

$\displaystyle\min_{\beta,z\in\mathcal{R}^{p},y\in\mathcal{R}^{n}}\Big{\\{}h(y)+\lambda\|z\|_{1}\
\Big{|}\ X\beta-y-b=0,\ \beta-z=0\Big{\\}}.$ (14)

Given $\sigma>0$, the augmented Lagrangian function corresponding to the
problem (14) takes the following form

$\displaystyle\mathcal{L}_{\sigma}(\beta,y,z;u,v)$ $\displaystyle:=$
$\displaystyle h(y)+\lambda\|z\|_{1}+\frac{\sigma}{2}\|X\beta-
y-b+\sigma^{-1}u\|^{2}-\frac{2}{\sigma}\|u\|^{2}$
$\displaystyle+\frac{\sigma}{2}\|\beta-z+\sigma^{-1}v\|-\frac{2}{\sigma}\|v\|^{2}.$

Therefore, we describe the details of the ADMM below.

Algorithm 4 (ADMM for the problem (1)):  Given
$\rho\in(0,\frac{1+\sqrt{5}}{2})$, $\sigma>0$, choose
$(y^{0},z^{0},u^{0},v^{0})\in\mathcal{R}^{n}\times\mathcal{R}^{p}\times\mathcal{R}^{n}\times\mathcal{R}^{p}$.
Set $k=0$ and iterate: Step 1. Compute $\displaystyle\beta^{k+1}$
$\displaystyle=$
$\displaystyle\mathop{\mbox{argmin}}_{\beta\in\mathcal{R}^{p}}\Big{\\{}\mathcal{L}_{\sigma}(\beta,y^{k},z^{k};u^{k},v^{k})\Big{\\}}$
$\displaystyle=$
$\displaystyle\left(I_{p}+X^{T}X\right)^{-1}\left(z^{k}-\sigma^{-1}v^{k}+X^{T}(y^{k}+b-\sigma^{-1}u^{k})\right);$
$\displaystyle(y^{k+1},z^{k+1})$ $\displaystyle=$
$\displaystyle\mathop{\mbox{argmin}}_{y\in\mathcal{R}^{n},z\in\mathcal{R}^{p}}\Big{\\{}\mathcal{L}_{\sigma}(\beta^{k+1},y,z;u^{k},v^{k})\Big{\\}}$
$\displaystyle=$
$\displaystyle\left(\mbox{Prox}_{\sigma^{-1}h}(X\beta^{k+1}-b+\sigma^{-1}u^{k}),\mbox{Prox}_{\sigma^{-1}\lambda\|\cdot\|_{1}}(\beta^{k+1}+\sigma^{-1}v^{k})\right).$
Step 2. Update the variables $u$ and $v$. $\displaystyle
u^{k+1}=u^{k}+\rho\sigma(X\beta^{k+1}-y^{k+1}-b),\quad
v^{k+1}=v^{k}+\rho\sigma(\beta^{k+1}-z^{k+1}).$ If the corresponding stopping
criteria is satisfied, then terminate; otherwise set $k:=k+1$ and goto Step 1.

####  5.1.2 Solving the problem (1) by linear programming

The problem (1) is polyhedral and can be solved via linear programming. By
introducing auxiliary variables, the problem (1) can be equivalently written
as the following form

$\displaystyle\min_{\beta,\xi,\zeta}\left\\{\frac{1}{n(n-1)}\sum_{1\leq
i<j\leq n}(\xi^{+}_{ij}+\xi^{-}_{ij})+\lambda\sum_{k=1}^{p}\zeta_{k}\right\\}$
$\displaystyle s.t.\
\xi^{+}_{ij}-\xi^{-}_{ij}=(b_{i}-X_{i}\beta)-(b_{j}-X_{j}\beta),\ 1\leq
i<j\leq n,$ $\displaystyle\quad\ \ \xi^{+}_{ij}\geq 0,\ \xi^{-}_{ij}\geq 0,\
1\leq i<j\leq n,$ $\displaystyle\quad\ \ \zeta_{k}\geq\beta_{k},\
\zeta_{k}\geq-\beta_{k},\ k=1,2\ldots,p.$

The above is a linear programming problem and can be solved by some existing
optimization software packages, e.g. Gurobi. In our comparison, we use the
Gurobi package to solve the above problem.

####  5.1.3 Numerical results for the tuning-free robust Lasso problem (1)

In this subsection, we show some performances of numerical results for the
tuning-free robust Lasso problem (1). In our comparison, we report the problem
(pbname), the number of samples ($n$), features ($p$) and nonzero elements
(nnz), lambda ($\lambda$), the relative KKT residual ($\eta_{kkt}$), the
primal objective value (pobj) and the running time (time) in the format of
“hours:minutes:seconds”. In the synthetic data sets, for a solution
$\widehat{\beta}$ obtained by an algorithm, we also report the $L_{1}$
estimation error $\|\widehat{\beta}-\dot{\beta}\|_{1}$ ($L_{1}$), the $L_{2}$
estimation error $\|\widehat{\beta}-\dot{\beta}\|_{2}$ ($L_{2}$), the model
error
$(\widehat{\beta}-\dot{\beta})^{T}\Sigma_{X}(\widehat{\beta}-\dot{\beta})$
(ME) with $\Sigma_{X}$ the covariance matrix of $X$, the number of false
positive variables (FP) which is the number of the noise covariates that are
selected in the model and the number of false negative variables (FN) which is
the number of the active variables that are not selected in the model. Except
the three columns of results including $\eta_{kkt}$, pobj and time, the other
columns of results are obtained by the PPMM algorithm.

Due to the excessive memory requirement of Gurobi, we first compare the PPMM
algorithm, ADMM and Gurobi on the synthetic datasets with $n=200$ and $p=800$.
The results are listed in Table 5.1.3.

Table 1: The performances of the PPMM, ADMM and Gurobi on the synthetic datasets ($n=200$, $p=1000$) for the tuning-free robust Lasso problem. In this table, “$a$”=PPMM, “$b$”=ADMM, “$c$”=Gurobi. pbname |  $\lambda$ |  nnz | $\eta_{kkt}$ | pobj | time |  $L_{1}$ |  $L_{2}$ |  ME |  FP |  FN
---|---|---|---|---|---|---|---|---|---|---
$n$; $p$ |  |  |  $a\ |\ b\ |\ c$ |  $a\ |\ b\ |\ c$ |  $a\ |\ b\ |\ c$ |  |  |  |  |
exmp1 |  0.271 |  4 |  2.3-7 $|$ 1.0-6 $|$ 2.2-10 |  1.4990+0 $|$ 1.4990+0 $|$ 1.4990+0 |  02 $|$ 1:27 $|$ 1:14 |  2.7-1 |  1.7-1 |  3.2-2 |  1 |  0
200;1000 |  |  |  |  |  |  |  |  |  |
exmp2 |  0.282 |  6 |  4.8-7 $|$ _2.0-6_ $|$ 3.3-10 |  1.9116+0 $|$ 1.9117+0 $|$ 1.9116+0 |  02 $|$ 1:54 $|$ 45 |  1.4+0 |  7.1-1 |  5.8-1 |  3 |  0
200;1000 |  |  |  |  |  |  |  |  |  |
exmp3 |  0.267 |  5 |  9.8-7 $|$ 1.0-6 $|$ 8.5-10 |  2.2355+0 $|$ 2.2355+0 $|$ 2.2355+0 |  02 $|$ 52 $|$ 53 |  2.0+0 |  1.1+0 |  2.0+0 |  2 |  0
200;1000 |  |  |  |  |  |  |  |  |  |
exmp4 |  0.272 |  4 |  3.6-7 $|$ _1.3-4_ $|$ 1.1-9 |  3.6425+0 $|$ 3.6437+0 $|$ 3.6425+0 |  02 $|$ 2:01 $|$ 53 |  4.6+0 |  2.5+0 |  1.2+1 |  1 |  0
200;1000 |  |  |  |  |  |  |  |  |  |
exmp5 |  0.276 |  7 |  6.3-7 $|$ _3.4-6_ $|$ 3.4-10 |  2.1721+0 $|$ 2.1721+0 $|$ 2.1721+0 |  02 $|$ 2:17 $|$ 1:04 |  1.9+0 |  1.0+0 |  1.9+0 |  4 |  0
200;1000 |  |  |  |  |  |  |  |  |  |
exmp6 |  0.270 |  9 |  9.9-7 $|$ _1.1-5_ $|$ 1.4-9 |  2.9175+0 $|$ 2.9177+0 $|$ 2.9175+0 |  02 $|$ 2:26 $|$ 1:29 |  2.5+0 |  1.2+0 |  2.0+0 |  6 |  0
200;1000 |  |  |  |  |  |  |  |  |  |
exmp7 |  0.266 |  94 |  9.9-7 $|$ _1.4-5_ $|$ 2.0-9 |  7.7347+0 $|$ 7.7352+0 $|$ 7.7347+0 |  06 $|$ 2:42 $|$ 1:09 |  1.4+1 |  2.0+0 |  6.9-1 |  72 |  3
200;1000 |  |  |  |  |  |  |  |  |  |
exmp8 |  0.267 |  88 |  9.5-7 $|$ _5.6-5_ $|$ 5.8-10 |  7.8796+0 $|$ 7.8816+0 $|$ 7.8796+0 |  05 $|$ 3:00 $|$ 53 |  1.9+1 |  2.9+0 |  2.3+0 |  68 |  5
200;1000 |  |  |  |  |  |  |  |  |  |
exmp9 |  0.268 |  79 |  4.5-7 $|$ _1.5-5_ $|$ 3.4-9 |  8.2800+0 $|$ 8.2804+0 $|$ 8.2800+0 |  05 $|$ 3:04 $|$ 1:01 |  2.4+1 |  3.5+0 |  4.4+0 |  59 |  5
200;1000 |  |  |  |  |  |  |  |  |  |
exmp10 |  0.269 |  59 |  5.3-7 $|$ _1.1-5_ $|$ 9.8-10 |  9.2758+0 $|$ 9.2760+0 $|$ 9.2758+0 |  04 $|$ 3:07 $|$ 50 |  3.9+1 |  6.3+0 |  2.4+1 |  49 |  15
200;1000 |  |  |  |  |  |  |  |  |  |
exmp11 |  0.272 |  67 |  4.0-7 $|$ _2.0-5_ $|$ 3.4-9 |  8.5739+0 $|$ 8.5747+0 $|$ 8.5739+0 |  04 $|$ 2:56 $|$ 1:01 |  1.9+1 |  3.0+0 |  3.5+0 |  45 |  3
200;1000 |  |  |  |  |  |  |  |  |  |
exmp12 |  0.259 |  58 |  7.8-7 $|$ _4.7-5_ $|$ 5.2-10 |  9.0375+0 $|$ 9.0383+0 $|$ 9.0374+0 |  04 $|$ 2:52 $|$ 52 |  4.0+1 |  6.2+0 |  1.7+1 |  50 |  17
200;1000 |  |  |  |  |  |  |  |  |  |

In the subsequent large-scale experiments, we list the performances of the
comparisons in Tables LABEL:table-compare-random-convex and LABEL:table-
compare-KEEL-convex. The symbol “–” in Tables LABEL:table-compare-random-
convex and LABEL:table-compare-KEEL-convex means that Gurobi fails to solve
the problem due to excessive memory requirement. In the implementation of the
ADMM, we apply the Sherman-Morrison-Woodbury formula [17] if it is necessary,
depending on the size of $n$ and $p$. We solve the linear system either by the
Cholesky factorization or by an iterative solver such as the preconditioned
conjugate gradient (PCG) method. From the comparisons, we can see that the
PPMM algorithm can obtain all the solutions with a desired accuracy
efficiently. When the sample size $n$ and the number of features $p$ are
small, Gurobi can also solve the problem with high accuracy but with much more
time. For high dimensional problems, Gurobi fails with too much memory
consuming.

###  5.2 Numerical experiments for the nonconvex tuning-free robust regression
problem

In this section, we implement some experiments for the tuning-free robust
regression problems with nonconvex regularizers. As we can see from Section
5.1.3 that Gurobi fails to solve high dimensional problems and it is time
consuming. Although we can reformulate each convex subproblem as a linear
programming problem and solve it by Gurobi, we only compare the the
performances of the ADMM and our PPMM algorithm for the nonconvex tuning-free
robust regression problem (2). We adopt the relative KKT residual

$\displaystyle\hat{\eta}_{kkt}:=\max\left\\{\frac{\|y-\mbox{Prox}_{h}(u+y)\|}{1+\|y\|},\frac{\|x-\mbox{Prox}_{\lambda
q_{1}-q_{2}}(x-A^{T}u)\|}{1+\|x\|},\frac{\|Ax-y-b\|}{1+\|y\|}\right\\}$

to measure the accuracy of the performance.

We first implement Stage 1 of our PPMM algorithm with $\eta_{kkt}<10^{-4}$ to
generate an initial point for Stage 2. All of the performances are stopped if
the relative KKT residual $\hat{\eta}_{kkt}<10^{-6}$. In addition, all of the
algorithms are also stopped when they reach the pre-set maximum number of
iterations (200 for Stage 2 of the PPMM algorithm and 40000 for the ADMM) or
the pre-set maximum running time of 4 hours. For each synthetic data set, we
generate 2000 observations for the training data set and 400 observations for
the validation data set. By applying the PPMM algorithm, we fit the model on
the training data set and use the validation data set to select the tuning
parameter $\lambda$. For each KEEL data set, we use the PPMM algorithm and
adopt a tenfold cross validation to find the parameter $\lambda$.

####  5.2.1 The ADMM for the nonconvex tuning-free robust regression problem

In this subsection, we describe the ADMM for solving the nonconvex tuning-free
robust regression problem. By introducing auxiliary variables, we can
reformulate the problem (2) as below

$\displaystyle\min_{\beta,z\in\mathcal{R}^{p},y\in\mathcal{R}^{n}}\Big{\\{}h(y)+\lambda
q_{1}(z)-q_{2}(z)\ \Big{|}\ X\beta-y-b=0,\ \beta-z=0\Big{\\}}.$ (15)

Given $\sigma>0$ and $(u,v)\in\mathcal{R}^{n}\times\mathcal{R}^{p}$, the
augmented Lagrangian function related to (15) is defined by

$\displaystyle\widehat{\mathcal{L}}_{\sigma}(\beta,y,z;u,v)$ $\displaystyle:=$
$\displaystyle h(y)+\lambda q_{1}(z)-q_{2}(z)+\frac{\sigma}{2}\|X\beta-
y-b+\sigma^{-1}u\|^{2}-\frac{1}{2\sigma}\|u\|^{2}$
$\displaystyle+\frac{\sigma}{2}\|\beta-z+\sigma^{-1}v\|^{2}-\frac{1}{2\sigma}\|v\|^{2}.$

We list the performance of the ADMM which is not guaranteed to converge due to
the nonconvexity in the following form.

Algorithm 5 (ADMM for the problem (2)):  Given $\sigma>0$, choose
$(y^{0},z^{0},u^{0},v^{0})\in\mathcal{R}^{n}\times\mathcal{R}^{p}\times\mathcal{R}^{n}\times\mathcal{R}^{p}$.
Set $k=0$ and iterate: Step 1. Compute $\displaystyle\beta^{k+1}$
$\displaystyle=$
$\displaystyle\mathop{\mbox{argmin}}_{\beta\in\mathcal{R}^{p}}\Big{\\{}\widehat{\mathcal{L}}_{\sigma}(\beta,y^{k},z^{k};u^{k},v^{k})\Big{\\}}$
$\displaystyle=$
$\displaystyle\left(I_{p}+X^{T}X\right)^{-1}\left(z^{k}-\sigma^{-1}v^{k}+X^{T}(y^{k}+b-\sigma^{-1}u^{k})\right);$
$\displaystyle(y^{k+1},z^{k+1})$ $\displaystyle=$
$\displaystyle\mathop{\mbox{argmin}}_{y\in\mathcal{R}^{n},z\in\mathcal{R}^{p}}\Big{\\{}\widehat{\mathcal{L}}_{\sigma}(\beta^{k+1},y,z;u^{k},v^{k})\Big{\\}}$
$\displaystyle=$
$\displaystyle\left(\mbox{Prox}_{\sigma^{-1}h}(X\beta^{k+1}-b+\sigma^{-1}u^{k}),\mbox{Prox}_{\sigma^{-1}(\lambda
q_{1}-q_{2})}(\beta^{k+1}+\sigma^{-1}v^{k})\right).$ Step 2. Update the
variables $u$ and $v$. $\displaystyle
u^{k+1}=u^{k}+\sigma(X\beta^{k+1}-y^{k+1}-b),\quad
v^{k+1}=v^{k}+\sigma(\beta^{k+1}-z^{k+1}).$ If the corresponding stopping
criterion is satisfied, then terminate; otherwise set $k:=k+1$ and return to
Step 1.

####  5.2.2 Numerical results

We first implement the numerical results for the tuning-free robust regression
problem with the SCAD regularization. The comparison results are listed in
Tables LABEL:table-compare-random-SCAD and LABEL:table-compare-SCAD-KEEL.
Comparing with the $\ell_{1}$ regularizer, the SCAD regularizer achieves a
sparse estimation with a fewer number of nonzeros for the synthetic type of
data set. The corresponding solution is much closer to $\dot{\beta}$ with less
values of $L_{1}$, $L_{2}$, ME and FP. For the first 8 synthetic data sets,
the SCAD regularizer can recover the true nonzero positions. The ADMM can
solve some problems with the desired accuracy, but it takes much more time
than the PPMM algorithm. From the comparison, we can see that the PPMM
algorithm can solve all the listed problems efficiently.

The second nonconvex problem we tested is the tuning-free robust regression
problem with the MCP regularization. The performances are listed in Tables
LABEL:table-compare-random-MCP and LABEL:table-compare-KEEL-MCP. For the
synthetic datasets, we can see that the MCP regularizer can recover the true
nonzero positions for the first 8 examples. Comparing with the $\ell_{1}$
regularizer, the MCP regularizer also achieves a sparse estimation with a
fewer number of nonzeros for the synthetic type of data set. The corresponding
solution is also much closer to $\dot{\beta}$ with less values of $L_{1}$,
$L_{2}$, ME and FP. The PPMM algorithm can solve all these problems with high
accuracy. In some results, the objective value obtained by the ADMM is smaller
than that of the PPMM, but the solutions derived by the PPMM algorithm have
fewer nonzeros. For example, the nnz for the ADMM are 7845, 3613, 11594, 169,
1224, 42 for the datasets example 12, baseball5, concrete7, dee10, friedman10
and mortgage5, respectively. The reason may be that Stage 1 of our PPMM
algorithm generates a relative better sparse initial point.

##  6 Conclusion

In this paper, we introduce the PPMM algorithm for solving nonconvex tuning-
free robust regression problems. We obtain the solution by solving a sequence
of convex majorization-minimization problems. The SSN based PPA is applied to
solve each convex subproblem. We make full use of the structure of the Clarke
generalized Jacobian and solve each subproblem efficiently. We also prove that
the PPMM algorithm converges to a d-stationary point and analyze the
convergence rate based on the KL property of the problem.

## References

  * [1] M. Ahn, J.S. Pang, and J. Xin, Difference-of-convex learning: directional stationarity, optimality, and sparsity, SIAM Journal on Optimization, 27 (2017), 1637–1665.
  * [2] J. Alcalá-Fdez, A. Fernandez, J. Luengo, J. Derrac, S. García, L. Sánchez, and F. Herrera, KEEL Data-Mining Software Tool: Data Set Repository, Integration of Algorithms and Experimental Analysis Framework, Journal of Multiple-Valued Logic and Soft Computing, 17 (2011), 255–287
  * [3] A. Beck, and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences, 2 (2009), 183–202.
  * [4] P.C. Bellec, G. Lecu$\acute{e}$, and A.B. Tsybakov, Slope meets Lasso: improved oracle bounds and optimality, The Annals of Statistics, 46 (2018), 3603–3642.
  * [5] A. Belloni, V. Chernozhukov, and L. Wang, Square-root Lasso: pivotal recovery of sparse signals via conic programming, Biometrika, 98 (2011), 791–806.
  * [6] M.J. Best, and N. Chakravarti, Active set algorithms for isotonic regression; a unifying framework, Mathematical Programming, 47 (1990), 425–439.
  * [7] J. Bolte, and E. Pauwels, Majorization-minimization procedures and convergence of SQP methods for semi-algebraic and tame programs, Mathematics of Operations Research, 41 (2016), 442–465.
  * [8] J. Bolte, S. Sabach, and M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Mathematical Programming, 146 (2014), 459–494.
  * [9] Y. Cui, J.S. Pang, and B. Sen, Composite difference-max programs for modern statistical estimation problems, SIAM Journal on Optimization, 28 (2018), 3344–3374.
  * [10] A. Derumigny, Improved bounds for square-root Lasso and square-root Slope, Electronic Journal of Statistics, 12 (2018), 741–766.
  * [11] J. Eckstein, and D.P. Bertsekas, On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators, Mathematical Programming, 55 (1992), 293–318.
  * [12] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, Least angle regression, The Annals of Statistics, 32 (2004), 407–499.
  * [13] F. Facchinei, and J.-S. Pang, Finite-dimensional Variational Inequalities and Com- plementarity Problems, Springer, New York, 2003.
  * [14] J.Q. Fan, and R.Z. Li, Variable selection via nonconcave penalized likelihood oracle property, Journal of the American Statistical Association, 96 (2001), 1348–1360.
  * [15] J.Q. Fan, and H. Peng, Nonconcave penalized likelihood with a diverging number of parameters, The Annals of Statistics, 32 (2004), 928–961.
  * [16] J. Flemming, Generalized Tikhonov regularization, basic theory and comprehensive results on convergence rates, PhD thesis, Fakultat fur Mathematik Technische Universitat Chemnitz, October 2011.
  * [17] G. Golub, and C.F. Van Loan, Matrix Computations, 3nd ed., Johns Hopkins University Press, Baltimore, MD, 1996.
  * [18] T.P. Hettmansperger, and J.W. McKean, Robust Nonparametric Statistical Methods, London: Arnold, 1998.
  * [19] L. Huang, J. Jia, B. Yu, B.G. Chun, P. Maniatis, and M. Naik, Predicting execution time of computer programs using sparse polynomial regression, in Advances in Neural Information Processing Systems, 2010, 883–891.
  * [20] J.-B. Hiriart-Urruty, J.-J. Strodiot, and V. H. Nguyen, Generalized Hessian matrix and second-order optimality conditions for problems with $C^{1,1}$ data, Applied Mathematics and Optimization, 11 (1984), 43–56.
  * [21] L.A. Jaeckel, Estimating regression coefficients by minimizing the dispersion of the residuals, The Annals of Mathematical Statistics, 43 (1972), 1449–1458.
  * [22] S.-J. Kim, K. Koh, M. Lustig, S. Boyd, and D. Gorinevsky, An interior-point method for large-scale $\ell_{1}$-regularized least squares, IEEE Journal of Selected Topics in Signal Processing, 1 (2007), 606–617.
  * [23] H.A. Le Thi, T. Pham Dinh, H.M. Le, and X.T. Vo, DC approximation approaches for sparse optimization, European Journal of Operations Research, 244 (2015), 26–46.
  * [24] X.D. Li, D.F. Sun, and K.-C. Toh, A highly efficient semismooth Newton augmented Lagrangian method for solving Lasso problems, SIAM Journal on Optimization, 28 (2018), 433–458.
  * [25] X.D. Li, D.F. Sun, and K.-C. Toh, On efficiently solving the subproblems of a level-set method for fused Lasso problems, SIAM Journal on Optimization, 28 (2018), 1842–1866.
  * [26] X.D. Li, D.F. Sun, and K.-C. Toh, An asymptotically superlinearly convergent semismooth Newton augmented Lagrangian method for linear programming, SIAM Journal on Optimization, 30 (2020), 2410–2440.
  * [27] X.G. Li, T. Zhao, X.M. Yuan, and H. Liu, The flare package for high dimensional linear regression and precision matrix estimation in R, Journal of Machine Learning Research, 16 (2015), 553–557.
  * [28] M.X. Lin, Y.-J. Liu, D.F. Sun, and K.-C. Toh, Efficient sparse semismooth Newton methods for the clustered Lasso problem, SIAM Journal on Optimization, 29 (2019), 2026–2052.
  * [29] B. Martinet, Régularisation d‘inéquations variationnelles par approximations successives, Revue Française Informatique Recherche Opérationnelle, 4 (1970), 154–158.
  * [30] F.W. Meng, D.F. Sun, and G.Y. Zhao, Semismoothness of solutions to generalized equations and the Moreau-Yosida regularization, Mathematical Programming, 104 (2005), 561–581.
  * [31] R. Mifflin, Semismooth and semiconvex functions in constrained optimization, SIAM Journal on Control and Optimization, 15 (1977), 959–972.
  * [32] L.Q. Qi, and J. Sun, A nonsmooth version of Newton’s method, Mathematical programming, 58 (1993), 353–367.
  * [33] S.M. Robinson, Some continuity properties of polyhedral multifunctions, Mathematical Programming at Oberwolfach, 14(1981), 206–214.
  * [34] R.T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, NJ, 1970.
  * [35] R.T. Rockafellar, Augmented Lagrangians and applications of the proximal point algorithm in convex programming, Mathematics of Operations Research, 1 (1976), 97–116.
  * [36] R.T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM Journal on Control and Optimization, 14 (1976), 877–898.
  * [37] R.T. Rockafellar, and R.J.-B. Wets, Variational Analysis, Springer, New York, 1998.
  * [38] P.P. Tang, C.J. Wang, D.F. Sun, and K.-C. Toh, A sparse semismooth Newton based proximal majorization-minimization algorithm for nonconvex square-root-loss regression problems, Journal of Machine Learning Research, 21 (2020), 1–38.
  * [39] R. Tibshirani, Regression shrinkage and selection via the Lasso, Journal of the Royal Statistical Society. Series B, 58 (1996), 267–288.
  * [40] L. Wang , B. Peng , J. Bradic , R.Z. Li, and Y.N. Wu, A tuning-free robust and efficient approach to high-dimensional regression, Journal of the American Statistical Association, 115 (2020), 1700–1714.
  * [41] C.H. Zhang, Nearly unbiased variable selection under minimax concave penalty, The Annals of Statistics, 38 (2010), 894–942.
  * [42] H. Zou, and R.Z. Li, One-step sparse estimates in nonconcave penalized likelihood models, The Annals of Statistics, 36 (2008), 1509–1533.