An Extended Newton-type Algorithm for $\ell_2$-Regularized Sparse Logistic Regression and Its Efficiency for Classifying Large-scale Datasets

TL;DR

This paper introduces an extended Newton-type algorithm for $ ext{l}_2$-regularized sparse logistic regression, demonstrating its efficiency and superior performance on large-scale datasets compared to existing solvers.

Contribution

It develops a novel Newton-based method for sparsity constrained logistic regression with proven convergence properties and improved computational efficiency.

Findings

Achieves global and quadratic convergence.

Outperforms seven state-of-the-art solvers in experiments.

Efficiently handles large-scale datasets.

Abstract

Sparse logistic regression, as an effective tool of classification, has been developed tremendously in recent two decades, from its origination the $\ell_1$-regularized version to the sparsity constrained models. This paper is carried out on the sparsity constrained logistic regression by the Newton method. We begin with establishing its first-order optimality condition associated with a $\tau$-stationary point. This point can be equivalently interpreted as a system of equations which is then efficiently solved by the Newton method. The method has a considerably low computational complexity and enjoys global and quadratic convergence properties. Numerical experiments on random and real data demonstrate its superior performance when against seven state-of-the-art solvers.

Figures

Tables

Table 1. Table 1 : Details of eight real datasets.

			Training	Testing
Data name	$n$	$p$	$m_1$	$m_2$
arcene	100	10,000	100	0
colon-cancer	62	2,000	62	0
news20.binary	19,996	1,355,191	19,996	0
newsgroup	11,314	777,811	11,314	0
duke breast-cancer	42	7,129	38	4
leukemia	72	7,129	38	34
gisette	7,000	5,000	6,000	1,000
rcv1.binary	40,242	47,236	20,242	20,000

Table 2. Table 2 : Benchmark Methods.

Models	(1.4) with convex ${\varphi }_{\nu }$	(1.4) with non-convex ${\varphi }_{\nu }$	(1.6)
first-order	SLEP	APG, GIST	GraSP
second-order	IRLS-LARS	--	NTGP, GPGN

Table 3. Table 3 : Sparsity levels ‖ 𝐳 ‖ 0 subscript norm 𝐳 0 \|\mathbf{z}\|_{0} of eight methods for Example 4.1 with p = 20000 , s = 0.1 ​ p formulae-sequence 𝑝 20000 𝑠 0.1 𝑝 p=20000,s=0.1p .

$n/p$	$0.1$	$0.2$	0.3	0.4	0.5	0.6	0.7
NSLR, GPGN	2000	2000	2000	2000	2000	2000	2000
GraSP, NTGP
LARS	500	500	500	500	500	500	500
GIST	4403	4309	5274	7832	7913	8614	8904
APG	6138	5857	5720	6170	5574	5048	5043
SLEP	2076	2534	2980	3235	3498	3596	3873

Table 4. Table 4 : Average results for Example 4.2 .

	$s=0.05p,n=0.2p$				$s=0.1p,n=0.2p$
	$\ell \left(𝐳\right)$	SER	Time	${\Vert 𝐳\Vert }_0$	$\ell \left(𝐳\right)$	SER	Time	${\Vert 𝐳\Vert }_0$
	$p=10000$
NSLR	3.2e-10	0.00e-0	0.436	500	1.1e-10	0.00e-0	0.918	1000
GPGN	1.09e-7	0.00e-0	4.192	500	6.84e-8	0.00e-0	6.962	1000
GraSP	1.20e-2	5.00e-3	17.10	500	7.41e-3	1.70e-3	17.01	1000
NTGP	4.37e-3	0.00e-0	33.54	500	2.18e-3	0.00e-0	13.53	1000
LARS	1.43e-1	1.25e-2	27.69	500	2.85e-1	4.90e-3	71.80	1000
GIST	4.81e-6	0.00e-0	12.54	2381	8.08e-6	0.00e-0	11.40	2974
APG	1.29e-4	0.00e-0	30.00	1407	1.25e-4	0.00e-0	27.82	5211
SLEP	1.75e-1	2.00e-3	6.560	811	1.32e-1	0.00e-0	10.84	1019
	$p=20000$
NSLR	1.6e-10	0.00e-0	2.522	1000	5.4e-11	0.00e-0	4.939	2000
GPGN	8.48e-8	0.00e-0	16.47	1000	5.91e-8	0.00e-0	15.83	2000
GraSP	1.29e-2	5.25e-3	62.22	1000	5.00e-3	1.50e-3	92.69	2000
NTGP	5.43e-3	0.00e-0	134.6	1000	2.30e-3	0.00e-0	54.50	2000
LARS	3.96e-1	5.43e-2	107.4	1000	4.21e-1	6.18e-2	117.1	1000
GIST	7.20e-7	0.00e-0	33.52	1542	9.98e-7	0.00e-0	54.55	4137
APG	5.06e-5	0.00e-0	24.94	1849	6.04e-5	0.00e-0	69.46	4323
SLEP	1.88e-1	3.25e-3	22.67	1511	1.45e-1	0.00e-0	34.59	2005
	$p=30000$
NSLR	1.1e-10	0.00e-0	7.364	1500	3.8e-11	0.00e-0	12.79	3000
GPGN	8.58e-8	0.00e-0	35.06	1500	4.09e-8	0.00e-0	51.49	3000
GraSP	2.28e-2	7.50e-3	181.9	1500	8.96e-3	2.83e-3	208.1	3000
NTGP	5.36e-3	0.00e-0	307.1	1500	2.32e-3	0.00e-0	125.1	3000
LARS	4.59e-1	1.03e-1	205.7	1000	4.99e-1	1.21e-1	206.7	1000
GIST	2.76e-6	0.00e-0	121.1	6124	7.97e-7	0.00e-0	150.9	6278
APG	9.40e-4	1.67e-4	206.7	7460	3.11e-4	0.00e-0	247.1	6450
SLEP	1.96e-1	2.33e-3	54.20	2392	1.49e-1	1.67e-4	80.19	3009

Table 5. Table 5 : Results for Example 4.3 with m 2 = 0 subscript 𝑚 2 0 m_{2}=0 .

	$\ell \left(𝐳\right)$	SER	Time	${\Vert 𝐳\Vert }_0$	$\ell \left(𝐳\right)$	SER	Time	${\Vert 𝐳\Vert }_0$
Data	Arcene				colon-cancer
NSLR	4.57e-7	0.00e-0	0.141	60	1.90e-8	0.00e-0	0.082	20
GPGN	1.84e-5	0.00e-0	0.583	60	3.35e-5	0.00e-0	0.171	20
GraSP	1.88e-3	0.00e-0	1.233	60	3.03e-1	1.61e-2	0.264	20
NTGP	1.09e-1	1.00e-2	2.812	60	4.51e-3	0.00e-0	0.192	20
LARS	7.98e-2	0.00e-0	0.974	60	2.84e-2	0.00e-0	0.231	30
GIST	3.99e-7	0.00e-0	5.644	134	6.22e-7	0.00e-0	0.133	41
APG	1.84e-5	0.00e-0	3.341	430	1.79e-7	0.00e-0	0.234	20
SLEP	1.05e-1	0.00e-0	3.663	66	1.68e-1	4.84e-2	0.142	23
Data	news20.binary				newsgroup
NSLR	1.11e-2	3.50e-3	3.236	2500	1.46e-2	8.00e-3	1.365	3000
GPGN	2.94e-2	8.00e-3	25.44	2500	5.17e-2	1.20e-2	30.71	3000
GraSP	2.46e-2	9.25e-3	205.7	2500	5.08e-1	4.75e-2	33.06	3000
NTGP	7.94e-2	6.55e-3	93.20	2500	3.11e-1	5.80e-2	13.07	3000
LARS	$--$	$--$	$--$	$--$	$--$	$--$	$--$	$--$
GIST	3.18e-2	6.60e-3	27.57	4091	6.84e-2	1.52e-2	10.55	3017
APG	4.18e-2	1.24e-2	37.15	3869	5.12e-2	1.77e-2	10.62	3257
SLEP	1.49e-1	3.34e-2	43.70	5299	2.60e-1	4.60e-2	23.89	5138

Table 6. Table 6 : Results for Example 4.3 with m 2 > 0 subscript 𝑚 2 0 m_{2}>0 .

	$\ell \left(𝐳\right)$	$\ell \left(𝐳\right)$-test	SER-train	SER-test	Time	${\Vert 𝐳\Vert }_0$
Data	duke breast-cancer
NSLR	4.10e-9	2.45e-6	0.00e-0	0.00e-0	0.113	100
GPGN	1.31e-5	2.55e-3	0.00e-0	0.00e-0	0.442	100
GraSP	3.57e-3	1.82e-3	0.00e-0	0.00e-0	0.502	100
NTGP	1.21e-5	1.20e-4	0.00e-0	0.00e-0	0.641	100
LARS	1.01e-4	1.21e-4	0.00e-0	0.00e-0	0.613	37
GIST	6.64e-9	1.60e-1	0.00e-0	0.00e-0	0.544	614
APG	6.35e-7	2.68e-7	0.00e-0	0.00e-0	0.742	136
SLEP	1.93e-3	1.04e-2	0.00e-0	0.00e-0	0.434	203
Data	leukemia
NSLR	3.09e-6	7.22e-2	0.00e-0	0.00e-0	0.113	150
GPGN	1.29e-5	2.71e-0	0.00e-0	1.47e-1	0.351	150
GraSP	3.40e-3	5.08e-1	0.00e-0	1.47e-1	0.562	150
NTGP	4.22e-4	2.11e-1	0.00e-0	5.88e-2	0.754	150
LARS	6.07e-4	1.11e-1	0.00e-0	8.82e-2	1.052	37
GIST	5.68e-3	1.82e-1	0.00e-0	1.18e-1	0.423	295
APG	1.05e-4	3.63e-1	0.00e-0	8.82e-2	0.542	1066
SLEP	1.71e-1	2.85e-1	0.00e-0	2.94e-2	0.581	269
Data	gisette
NSLR	2.88e-6	6.51e-1	0.00e-0	4.30e-2	1.221	500
GPGN	2.25e-4	2.63e-1	0.00e-0	4.60e-2	1.412	500
GraSP	1.51e-4	2.30e-1	0.00e-0	4.82e-2	2.433	500
NTGP	1.17e-3	9.58e-1	0.00e-0	4.18e-2	2.554	500
LARS	1.63e-1	1.39e-0	1.00e-3	4.18e-2	8.262	500
GIST	2.40e-4	1.02e-0	0.00e-0	4.30e-2	1.783	1303
APG	3.72e-4	1.31e-0	0.00e-0	4.90e-2	1.641	907
SLEP	1.92e-2	8.30e-1	0.00e-0	4.97e-2	2.042	1569
Data	rcv1.binary
NSLR	5.82e-2	2.01e-1	1.95e-2	5.45e-2	3.712	1000
GPGN	2.90e-2	2.35e-1	8.05e-3	5.58e-2	6.811	1000
GraSP	3.09e-1	1.98e-0	4.15e-2	9.51e-2	21.72	1000
NTGP	7.48e-2	1.37e-1	1.06e-2	4.64e-2	4.471	1000
LARS	2.27e-1	2.61e-1	5.13e-1	5.46e-1	33.53	1000
GIST	3.44e-2	1.41e-1	8.20e-3	4.85e-2	4.571	1545
APG	2.51e-2	2.78e-1	7.31e-3	6.00e-2	7.571	1537
SLEP	1.24e-1	1.65e-1	3.02e-2	5.23e-2	8.453	3527

Equations

$$\mathrm{Pr}\left(y|\mathbf{x},\mathbf{z}\right)={\left(1+e^{-⟨\mathbf{x},\mathbf{z}⟩}\right)}^{-1},$$

$$\ell (\mathbf{z}):=\frac{1}{n}\sum _{i=1}^n\left(\ln \left(1+e^{⟨{\mathbf{x}}_i,\mathbf{z}⟩}\right)-y_i⟨{\mathbf{x}}_i,\mathbf{z}⟩\right).$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^p}{\min }\ell (\mathbf{z})+\nu \Vert \mathbf{z}{\Vert }_1,$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^p}{\min }{\ell }_{\varphi }(\mathbf{z}):=\ell (\mathbf{z})+{\varphi }_{\nu }(\mathbf{z}),$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^p}{\min }\ell (\mathbf{z}),s.t.\Vert \mathbf{z}{\Vert }_0\le s,$$

$$\min f(\mathbf{z}):=\ell (\mathbf{z})+(\lambda /2)\Vert \mathbf{z}{\Vert }_2^2,s.t.\Vert \mathbf{z}{\Vert }_0\le s,$$

$${\lambda }_{\mathrm{x}}:={\lambda }_{\max }(X^{T}X)/(4n).$$

$$\displaystyle\begin{array}[]{lll}f(\mathbf{z})&\geq&f(\mathbf{z}^{\prime})+\langle\nabla f(\mathbf{z}^{\prime}),\mathbf{z}-\mathbf{z}^{\prime}\rangle+({\lambda}/{2})\|\mathbf{z}-\mathbf{z}^{\prime}\|^{2},\\ f(\mathbf{z})&\leq&f(\mathbf{z}^{\prime})+\langle\nabla f(\mathbf{z}^{\prime}),\mathbf{z}-\mathbf{z}^{\prime}\rangle+({L}/{2})\|\mathbf{z}-\mathbf{z}^{\prime}\|^{2}.\end{array}$$

$$\Vert \nabla f(\mathbf{z})-\nabla f({\mathbf{z}}^{\prime })\Vert \le L\Vert \mathbf{z}-{\mathbf{z}}^{\prime }{\Vert }^2,$$

$$\Vert {\nabla }^{2}f(\mathbf{z})-{\nabla }^{2}f({\mathbf{z}}^{\prime })\Vert \le M\Vert \mathbf{z}-{\mathbf{z}}^{\prime }{\Vert }^2.$$

$${\Pi }_s(\mathbf{z}):={\mathrm{argmin}}_{\mathbf{x}\in {\mathbb{R}}^p}\{\Vert \mathbf{z}-\mathbf{x}\Vert :\Vert \mathbf{x}{\Vert }_0\le s\},$$

$$\mathbf{z}\in {\Pi }_s\left(\mathbf{z}-\tau \nabla f(\mathbf{z})\right).$$

$$\displaystyle\|\mathbf{z}\|_{0}\leq s,~{}~{}~{}~{}\tau\left|(\nabla f(\mathbf{z}))_{i}\right|\left\{\begin{array}[]{lll}=&0,&i\in\mathrm{supp}(\mathbf{z}),\\ \leq&[\mathbf{z}]^{\downarrow}_{s},&i\notin\mathrm{supp}(\mathbf{z}).\end{array}\right.$$

$$\displaystyle\tau\left|(\nabla f(\mathbf{z}^{*}))_{i}\right|\left\{\begin{array}[]{lll}=&0,&i\in\mathrm{supp}(\mathbf{z}^{*})=:\alpha_{*},\\ \leq&[\mathbf{z}^{*}]^{\downarrow}_{s},&i\notin\mathrm{supp}(\mathbf{z}^{*}).\end{array}\right.$$

$$\displaystyle\begin{array}[]{lll}2f(\mathbf{z})-2f(\mathbf{z}^{*})&\geq&2\langle\nabla f(\mathbf{z}^{*}),\mathbf{z}-\mathbf{z}^{*}\rangle+{\lambda}\|\mathbf{z}-\mathbf{z}^{*}\|^{2}\\ &=&2\left\langle(\nabla f(\mathbf{z}^{*}))_{\overline{\alpha}_{*}},(\mathbf{z}-\mathbf{z}^{*})_{\overline{\alpha}_{*}}\right\rangle+{\lambda}\|\mathbf{z}-\mathbf{z}^{*}\|^{2}\\ &\overset{\eqref{agradient-*}}{=}&{\lambda}\|\mathbf{z}-\mathbf{z}^{*}\|^{2}.\end{array}$$

$$\Vert {\mathbf{z}}^{\ast }-({\mathbf{z}}^{\ast }-\tau \nabla f({\mathbf{z}}^{\ast })){\Vert }^2\le \Vert \mathbf{z}-({\mathbf{z}}^{\ast }-\tau \nabla f({\mathbf{z}}^{\ast })){\Vert }^2.$$

$$2⟨\nabla f({\mathbf{z}}^{\ast }),\mathbf{z}-{\mathbf{z}}^{\ast }⟩\ge -(1/\tau )\Vert \mathbf{z}-{\mathbf{z}}^{\ast }{\Vert }^2.$$

$$2f(\mathbf{z})\ge 2f({\mathbf{z}}^{\ast })+(\lambda -1/\tau )\Vert \mathbf{z}-{\mathbf{z}}^{\ast }{\Vert }^2,$$

$${\Sigma }_s(\mathbf{z}):=\{\alpha \in [p]:|\alpha |=s,|z_i|\ge |z_j|,\forall i\in \alpha ,j\in \overline{\alpha }\}.$$

$$\displaystyle\Pi_{s}(\mathbf{z})=\Big{\{}\left(\mathbf{z}_{\alpha}^{T}~{}0\right)^{T}:~{}\alpha\in\Sigma_{s}(\mathbf{z})\Big{\}}.$$

$$(\nabla f(\mathbf{z}){)}_{\alpha }=0,{\mathbf{z}}_{\overline{\alpha }}=0.$$

$$H^k:={\nabla }^{2}f({\mathbf{z}}^k),{\mathbf{g}}^k:=\nabla f({\mathbf{z}}^k).$$

$$\displaystyle\begin{pmatrix}H^{k}_{\alpha\alpha}&H^{k}_{\alpha\overline{\alpha}}\\ 0&I_{p-s}\end{pmatrix}\mathbf{d}^{k}=-\left(\begin{array}[]{c}{\boldsymbol{g}}^{k}_{\alpha}\\ \mathbf{z}^{k}_{\overline{\alpha}}\end{array}\right)=:-\theta_{\alpha}^{k},$$

$$\{\begin{array}{ll}{H}_{\alpha \alpha }^k{\mathbf{d}}_{\alpha }^k& =H_{\alpha \overline{\alpha }}^k{\mathbf{z}}_{\overline{\alpha }}^k-{\mathbf{g}}_{\alpha }^k,\\ {\mathbf{d}}_{\overline{\alpha }}^k& =-{\mathbf{z}}_{\overline{\alpha }}^k.\end{array}$$

$$\Vert {\mathbf{z}}^{k+1}{\Vert }_0=\Vert {\mathbf{z}}_{\alpha }^{k+1}{\Vert }_0\le |\alpha |=s.$$

$${\mathbf{z}}^k(\sigma ):=\left(\begin{array}{c}{\mathbf{z}}_{\alpha }^k+\sigma {\mathbf{d}}_{\alpha }^k\\ {\mathbf{z}}_{\overline{\alpha }}^k+{\mathbf{d}}_{\overline{\alpha }}^k\end{array}\right)=\left(\begin{array}{c}{\mathbf{z}}_{\alpha }^k+\sigma {\mathbf{d}}_{\alpha }^k\\ 0\end{array}\right).$$

$$2f({\mathbf{z}}^k(c^r))\le 2f({\mathbf{z}}^k)+c^r⟨{\mathbf{g}}^k,{\mathbf{d}}^k⟩.$$

$$\beta \in {\Sigma }_s({\mathbf{z}}^{k-1}-\tau {\mathbf{g}}^{k-1}).$$

$$\mathrm{supp}({\mathbf{z}}^k)\subseteq \beta ,\mathrm{supp}({\mathbf{z}}^{k+1})\subseteq \alpha .$$

$$\displaystyle-\mathbf{d}^{k}_{\overline{\alpha}}=\mathbf{z}^{k}_{\overline{\alpha}}=\left(\begin{array}[]{c}\mathbf{z}^{k}_{\beta\cap\overline{\alpha}}\\ 0\end{array}\right)=\left(\begin{array}[]{c}\mathbf{z}^{k}_{\beta\setminus{\alpha}}\\ 0\end{array}\right)=\left(\begin{array}[]{c}\mathbf{z}^{k}_{\gamma}\\ 0\end{array}\right).$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

As one of effective tools of classification, logistic regression has its high reputation with extensive applications ranging from machine learning, data mining, pattern recognition, medical science to statistics. It describes the relationship between a sample data $\mathbf{x}$ and its associated binary response/label $y\in\{0,1\}$ through the conditional probability

[TABLE]

where ${\rm Pr}(y|\mathbf{x},\mathbf{z})$ is the conditional probability of the label $y$, given the sample $\mathbf{x}$ and a parameter vector $\mathbf{z}$, and $\langle\mathbf{x},\mathbf{z}\rangle$ is the vector inner product. To find the maximum likelihood estimate of the parameter $\mathbf{z}$, a set of $n$ i.i.d. (independently and identically distributed) samples $\{(\mathbf{x}_{i},y_{i}),i=1,2,\ldots,n\}$ are first drawn, where $\mathbf{x}_{i}\in\mathbb{R}^{p}$ and $y_{i}\in\{0,1\}$, yielding a joint likelihood of the interested parameter/classifier $\mathbf{z}$. Then the maximum likelihood estimate is obtained by minimizing the classical logistic regression loss function,

[TABLE]

The logistic loss function is strictly convex and thus admits a unique minimizer provided that the sample matrix is full row rank. Therefore, the minimization performs relatively well when the number of samples is larger than the number of features, i.e., $n\geq p$. But, the case $n<p$ may lead to an over-fitting: the solved classifier through minimizing (1.2) well fits the model (making the loss sufficiently small) on training data but behaves poorly on unseen data.

On the one hand, the case $n<p$ occurs often in many real applications. For instance, one piece of gene expression data sample is made of thousands of genes whilst common medical equipments are only able to obtain very limited samples. In image processing, an image consists of large amounts of pixels, which is far more than the number of observed images. One the other hand, despite numerous features in those data, there is only a small portion that is of importance. For example, apart from the classification task, the micro-array data experiments also attempt to identify a small set of informative genes (to distinguish the tumour and the normal tissues) in each gene expression data so as to remove the irrelevant genes to simplify the inference. This naturally gives rise to the topic of the sparse logistic regression.

1.1 Sparse logistic regression

Sparse logistic regression (SLR) was originated from the $\ell_{1}$-regularized logistic regression in [1],

[TABLE]

where $\|\mathbf{z}\|_{1}$ is the $\ell_{1}$-norm and $\nu>0$. Under the help of $\ell_{1}$-regularization, this model is capable of rendering a sparse solution allowing for capturing key features among others. A vector is called sparse if only a few entries are non-zero and the rest are zeros. With the advance in sparse optimization in recent decade, (1.3) has been extensively extended to the following general model,

[TABLE]

where the regularized function $\phi_{\nu}(\mathbf{z}):\mathbb{R}^{p}\rightarrow\mathbb{R}$ is designed to pursue a sparse solution and associated with some given non-negative parameters $\nu$.

An alternative is to consider logistic regression with a sparsity constraint, which was first studied in [2, 3] separately and then well investigated in [4]. They perform the following sparsity constrained logistic regression

[TABLE]

where $\|\mathbf{z}\|_{0}$ is the $\ell_{0}$ pseudo norm of $\mathbf{z}$, counting the number of non-zero elements of $\mathbf{z}$. The discreteness of the sparsity constraint makes tackling this model NP-hard. Nevertheless, compared with the regularized model, the sparsity constrained version enjoys various appealing features, such as being penalty parameter-free, ease of sparsity controlling, and low computational complexity in terms of numerical computation and so forth.

Therefore, the generalization of the problem (1.5), where $\ell(\mathbf{z})$ is replaced by a more general function, has been thoroughly investigated in [5, 6] since it was first introduced by [2] and [7]. Particularly, in statistics, the model with the logistic loss function being replaced by the least squares of linear regression is the so-called best subspace/feature selection [8, 9, 10, 11, 12]. Those research bring fruitful results and provide a series of effective numerical tools to conquer the NP-hardness.

However, as stated in [2] that ‘one can achieve arbitrarily small loss values by tending the parameters to infinity along certain directions’ for (1.5), authors [2] suggests to address the following regularized model

[TABLE]

where $\lambda>0$ is a given penalty parameter. Now the objective function $f$ is strongly convex and thus (1.6) admits finitely many (local or global) bounded minimizers. So the work in this paper is carried out along with this model.

1.2 Methods of solving SLR

Since there is a vast body of methods that have been proposed to deal with the sparse optimization problems containing the SLR as a special case, we present a brief overview of methods that process the problems (1.4)-(1.6) directly.

Regularization methods

Most versions of the model (1.4) are unconstrained and continuous. Then generic optimization methods, known as the relaxation (regularization) methods from the perspective of optimization, are tractable. Dependent on the convexity of the penalty functions $\phi_{\nu}$, those methods can be summarized into two categories.

Convex regularizations are mainly associated with the usage of $\ell_{1}$-norm:

$\phi_{\nu}(\mathbf{z})=\nu\|\mathbf{z}\|_{1}$. Some earliest work can be traced back to [13, 14], where expectation maximization methods were developed. Later relevant work can be found in [15, 16, 17, 18, 19].

• 2.

$\phi_{\boldsymbol{\nu}}(\mathbf{z})=\nu_{1}\|\mathbf{z}\|^{2}_{2}+\nu_{2}\|\mathbf{z}\|_{1}$, where $\boldsymbol{\nu}:=(\nu_{1},\nu_{2})>0$. For this penalty, two powerful packages SLEP [20] and GLMNET [21, 22] have been created.

• 3.

$\phi_{\nu}(\mathbf{z})=\nu\|\mathbf{z}\|^{2}+\delta_{\|\mathbf{z}\|_{1}\leq t}(\mathbf{z})$, where $t>0$ is a given parameter, and $\delta_{\|\mathbf{z}\|_{1}\leq t}(\mathbf{z})=0$ if $\|\mathbf{z}\|_{1}\leq t$ and $+\infty$ otherwise. Such a problem can be addressed by Lassplore [23] or SLEP [20]. When $\nu=0$, the above model is the $\ell_{1}$ constrained logistic regression, which was addressed by IRLS-LARS in [24]. Here, LARS was adopt from [25].

Nonconvex regularizations differ slightly. In the early stage, scholars from statistics have proposed a number of excellent methods including the smoothly clipped absolute deviation (SCAD [26]), one step local linear approximation [27] and the group bridge method for multiple regression problems [28]. Then, a general iterative shrinkage and thresholding algorithm (GIST) has been proposed in [29]. Recently, the accelerated proximal gradient method (APG) in [30], the efficient hybrid optimization algorithm for non-convex regularized problems (HONOR) in [31] and the proximal Newton method based on the scheme of the difference of two convex functions in [32] are worth exploring.

Greedy Methods

An impressive body of approaches have been developed to solve the sparsity constrained models (1.5) or (1.6). The first work in [2] generalized the compressive sampling matching pursuit [33] to derive the gradient support pursuit (GraSP). Then authors in [34] adopted the orthogonal matching pursuit (OMP [35]) to develop a group OMP method. Other relating methods can be seen those in [36, 37, 38]. Very lately, three effective Newton type methods have been designed. They are the Newton greedy pursuit method NTGP in [39], greedy projected gradient-Newton method (GPGN [4]) and the fast Newton hard thresholding pursuit [40]. In particular, we would like to mention the methods, the zero-CW search method and the full-CW search method, proposed in [5]. Both methods first carefully search an index set $T$ and then solve a subproblem where the variable has support within $T$ to update the next point.

1.3 Our contributions

Those aforementioned methods have been testified to have the excellent numerical performance to deal with (1.5) or (1.6). However, only a very few of them established strong theoretical guarantees (such as global convergence property or quadratic convergence rate) from the perspective of deterministic optimization. Therefore, in this paper, we aim to develop a second-order method that possesses such strong theoretical guarantees. The main contributions are summarized as follows.

C1) We start with establishing the optimality condition of the model (1.6) by introducing a $\tau$-stationary point (see Definition 2.2 for more details) which turns out to be at least a locally optimal solution by Theorem 2.3. More importantly, a $\tau$-stationary point draws forth a system of equations (2.17) that makes the classic Newton method applicable.

C2) Differing with any of the above mentioned algorithms, we perform the Newton method on solving a system of equations (2.17), one of the optimality conditions of the problem (1.6). The proposed Newton method for SLR (NSLR for short) has a simple framework (see Algorithm 1) that makes its implementation easy and a low computational complexity per each iteration. Such a low computational complexity is due to a small-scale linear equation system with $s$ variables and $s$ equations being solved to update the Newton direction.

C3) It is worth mentioning that the standard Newton-type methods derive the directions for a fixed system of equations. However, in each iteration, the system of equations (2.17) varies when the index set $\alpha$ changes. Consequently, NSLR  updates Newton directions on unfixed systems of equations. Because of this, some common approaches to establish the convergence results of Newton-type methods for solving a fixed system of equations fail to be employed for NSLR. Nevertheless, we still show that the whole sequence generated by NSLR converges to a $\tau$-stationary point, at least a locally optimal solution. Moreover, the convergence enjoys a quadratic rate, well testifying the proposed method would perform extraordinarily theoretically.

C4) Finally, the efficiency of NSLR is demonstrated against seven state-of-the-art methods by solving a number of randomly generated and real datasets. The fitting accuracy and computational speed are very competitive. Especially, in high dimensional data setting, NSLR outperforms the others in terms of the computational time.

We note that there are some methods that also have a close link to $\tau$-stationary point, such as two methods in [5] and GPGN [4]. We would like to highlight the difference between them and NSLR. For the methods in [5], since an optimal solution to a subproblem needs to be found to update the next point in each step, the methods can terminate within finitely many steps. However, NSLR updates the next point by Newton direction with a line search scheme and has been shown to enjoy the global and quadratic convergence properties. For GPGN, the procedure IHT from [7] to update the next point dominates most steps, and Newton steps are imposed only when two consecutive points have the same support sets. It is shown to converge quadratically only when the solution has $s$ nonzeros but sublinearly otherwise. By contrast, NSLR always performs Newton step to update the next point and converges quadratically without additional assumptions. Moreover, differing from papers [5] and [4] where comprehensive optimality conditions have been investigated, the primary aim of this paper is to develop a Newton-type method and establish its convergence properties.

1.4 Organization and notation

This paper is organized as follows. To explore the optimality conditions of (1.6), the next section introduces the $\tau$-stationary point by Definition 2.2 and establishes its relationships with a local/global minimizer in Theorem 2.3. This $\tau$-stationary point is then equivalently transferred to a system of equations (2.17). Section 3 develops the method NSLR, an abbreviation for Newton method for SLR, which turns out to have a simple algorithmic framework and low computational complexity. The global and quadratic convergence properties of the method are then established. In Section 4, the superior performance of NSLR is demonstrated against some of the state-of-the-art solvers on randomly generated and real datasets in high dimensional scenarios. Concluding remarks are made in the last section.

We end this section by defining some notation employed throughout this paper. Let $X:=(\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n})^{T}\in\mathbb{R}^{n\times p}$ be the sample matrix and $\mathbf{y}=(y_{1},y_{2},\ldots,y_{n})^{T}\in\mathbb{R}^{n}$ be the response vector. For an index set $\alpha\subseteq[p]:=\{1,2,\ldots,p\}$, let $|\alpha|$ be the cardinality of $\alpha$ and $\overline{\alpha}:=[p]\setminus\alpha$ be the complementary set of $\alpha$. The support set of a vector $\mathbf{z}$ is denoted by $\mathrm{supp}(\mathbf{z}):=\{i\in[p]:z_{i}\neq 0\}$. We denote $[\mathbf{z}]^{\downarrow}_{i}$ the $i$th largest (in absolute) elements of $\mathbf{z}$. Write $\mathbf{z}_{\alpha}\in\mathbb{R}^{|\alpha|}$ as the sub-vector of $\mathbf{z}$ containing elements indexed on $\alpha$. Similarly, for a matrix $A\in\mathbb{R}^{p\times p}$, $A_{\alpha\beta}$ is the sub-matrix containing rows indexed on ${\alpha}$ and columns indexed on $\beta$, particularly, $A_{\alpha:}=A_{\alpha[p]}$ if $\beta=[p]$. Let $\|\cdot\|$ denote the Spectral norm for a matrix and Euclidean norm for a vector respectively. Furthermore, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ are the minimal and maximal eigenvalues of $A$.

2 Optimality

This section is devoted to investigate the optimality conditions of (1.6), before which we summarize some properties of $\ell(\mathbf{z})$ from [4].

Proposition 2.1** (Lemma 2.2-2.4, Lemma A.3 [4])**

The function $\ell(\mathbf{z})$ is twice continuously differentiable and has the following basic properties:

• i)

It is non-negative, convex and strongly smooth on $\mathbb{R}^{p}$ with a parameter

[TABLE]

• ii)

The gradient is Lipschitz continuous with the Lipschitz constant $\lambda_{\mathrm{x}}$.

• iii)

The Hessian matrix is Lipschitz continuous with the Lipschitz constant $M:=12\lambda_{\mathrm{x}}{\max}_{i\in[n]}\|\mathbf{x}_{i}\|_{1}$.

These properties of $\ell(\mathbf{z})$ are also enjoyed by the function $f(\mathbf{z})=\ell(\mathbf{z})+(\lambda/2)\|\mathbf{z}\|^{2}$. Since the proofs are easy, we only summarize them here. The function $f$ is strongly convex with a constant $\lambda$, and strongly smooth with a parameter $L:=\lambda+\lambda_{\mathrm{x}}$, namely, for any $\mathbf{z},\mathbf{z}^{\prime}\in\mathbb{R}^{p}$,

[TABLE]

The gradient is Lipschitz continuous with the Lipschitz constant $L$, namely,

[TABLE]

for any $\mathbf{z},\mathbf{z}^{\prime}\in\mathbb{R}^{p}$. The Hessian matrix is Lipschitz continuous with the Lipschitz constant $M$, namely, for any $\mathbf{z},\mathbf{z}^{\prime}\in\mathbb{R}^{p}$,

[TABLE]

When it comes to characterize the solutions of the problem $(\ref{SLR-L2})$, we need the projection of a vector $\mathbf{z}\in\mathbb{R}^{p}$ onto the feasible region defined by

[TABLE]

which sets all but $s$ largest absolute value components of $\mathbf{z}$ to zero. Since the right hand side may have multiple solutions, $\Pi_{s}(\mathbf{z})$ is a set. Based on the projection, we introduce the concept of the $\tau$-stationary point which is also known as the $L$ stationary point in [7, Definition 2.3].

Definition 2.2

[7, Definition 2.3]** A point $\mathbf{z}$ is called a $\tau$-stationary point of the problem $\eqref{SLR-L2}$ if there is a $\tau>0$ satisfying

[TABLE]

By [7, Lemma 2.2 ], $\mathbf{z}$ is a $\tau$-stationary point if and only if

[TABLE]

Based on the definition of the $\tau$-stationary point, our first main result is establishing its relationships with a locally/globally optimal solution to (1.6).

Theorem 2.3

The following results hold for the problem $(\ref{SLR-L2})$.

• i)

A global minimizer is a $\tau$-stationary point $\mathbf{z}^{*}$ for any $0<\tau<1/L.$

• ii)

A $\tau$-stationary point $\mathbf{z}^{*}$ for some $\tau>0$ is a unique local minimizer if $\|\mathbf{z}^{*}\|_{0}=s$ and a unique global minimizer if $\|\mathbf{z}^{*}\|_{0}<s$.

• iii)

A $\tau$-stationary point for some $\tau>1/\lambda$ is a unique global minimizer.

Proof i)  The proof is the same as that in [7, Theorem 2.2].

ii)  Let $\mathbf{z}^{*}$ be a $\tau$-stationary point for some $\tau>0$. Then we have (2.8), i.e.,

[TABLE]

For the case $\|\mathbf{z}^{*}\|_{0}=s$, consider a local region $N(\mathbf{z}^{*}):=\{\mathbf{z}:\|\mathbf{z}-\mathbf{z}^{*}\|<[\mathbf{z}^{*}]^{\downarrow}_{s}\}$. Then for any feasible point $\mathbf{z}\in N(\mathbf{z}^{*})$ and any $i\in\alpha_{*}$, we have $|z_{i}|\geq|z^{*}_{i}|-|z^{*}_{i}-z_{i}|>|z^{*}_{i}|-[\mathbf{z}^{*}]^{\downarrow}_{s}\geq 0,$ which means $\alpha_{*}\subseteq\mathrm{supp}(\mathbf{z})$. Since $\|\mathbf{z}\|_{0}\leq s=|\alpha_{*}|$, it holds $\alpha_{*}=\mathrm{supp}(\mathbf{z})$ for any feasible point $\mathbf{z}\in N(\mathbf{z}^{*})$, namely $\mathbf{z}_{\overline{\alpha}_{*}}=\mathbf{z}^{*}_{\overline{\alpha}_{*}}=0$. Then the strong convexity of $f$ in (2.3) leads to

[TABLE]

Thus $\mathbf{z}^{*}$ is a unique local minimizer of (1.6).

For the case $\|\mathbf{z}^{*}\|_{0}<s$, the condition (2.11) implies $\nabla f(\mathbf{z}^{*})=0$ due to $[\mathbf{z}^{*}]^{\downarrow}_{s}=0$. Then (2.15) is true for any $\|\mathbf{z}\|_{0}\leq s$. So $\mathbf{z}^{*}$ is a unique global minimizer of the problem (1.6).

iii)   Let $\mathbf{z}^{*}$ be a $\tau$-stationary point for some $\tau>0$. Then (2.5) and the definition of the projection $\Pi_{s}$ imply that

[TABLE]

This leads to

[TABLE]

The above condition together with the inequality in (2.15) derives

[TABLE]

which shows the unique global optimality of $\mathbf{z}^{*}$ if $\tau>1/\lambda$. \qed

We note that the necessary optimality condition in Theorem 2.3 i) is directly adopted from [7, Theorem 2.2] or [5, Theorem 5.3] where $\mathbb{B}=\mathbb{R}^{p}$. However, we also establish the sufficient optimality conditions, see Theorem 2.3 ii) and iii). The above relationships show that a $\tau$-stationary point is at least a unique locally optimal solution to the problem (1.6). This allows us to focus on a $\tau$-stationary point itself to pursue a ‘good’ solution. Therefore, we define a set

[TABLE]

Each element $\alpha$ in $\Sigma_{s}(\mathbf{z})$ coincides the indices of the first $s$ largest (in absolute) components of $\mathbf{z}$. Note that $\Sigma_{s}(\mathbf{z})$ may have multiple elements. For instance, $\mathbf{z}=(3,-2,2,1,0)^{T}$, $\Sigma_{3}(\mathbf{z})=\{\{1,2,3\}\}$ and $\Sigma_{2}(\mathbf{z})=\{\{1,2\},\{1,3\}\}.$ The notation allows us to rewrite $\Pi_{s}(\mathbf{z})$ as follows

[TABLE]

Then a point satisfying (2.5) can be interpreted as that there is an $\alpha\in\Sigma_{s}(\mathbf{z}-\tau\nabla f(\mathbf{z}))$ satisfying $\mathbf{z}_{\alpha}=(\mathbf{z}-\tau\nabla f(\mathbf{z}))_{\alpha}$ and $\mathbf{z}_{\overline{\alpha}}=0$, which is equivalent to

[TABLE]

Therefore, to find a $\tau$-stationary point of (1.6), one can seek for a solution to the equation system (2.17). This is summarized into the following theorem.

Theorem 2.4

A point $\mathbf{z}$ is a $\tau$-stationary point of (1.6) if and only if there is an $\alpha\in\Sigma_{s}(\mathbf{z}-\tau\nabla f(\mathbf{z}))$ satisfying (2.17).

3 Convergence Analysis

In this section, we turn our attention to solve the equations (2.17) to pursue a $\tau$-stationary point of the problem (1.6), at least a unique local minimizer.

Given a point $\mathbf{z}^{k}$, for notational convenience, let

[TABLE]

3.1 The framework

Suppose we have a point $\mathbf{z}^{k}$ computed already. Then we can pick an index set $\alpha$ from $\Sigma_{s}(\mathbf{z}^{k}-\tau{{\boldsymbol{g}}}^{k})$. For such a fixed index set $\alpha$, we apply Newton step on the equations (2.17) just once to derive the Newton direction by

[TABLE]

where $I_{p-s}$ is the $(p-s)$ order identity matrix. One can calculate $\mathbf{d}^{k}$ by

[TABLE]

Since $f$ is strongly convex, $H^{k}$ is non-singular and so are its any principal sub-matrices, i.e., $H^{k}_{\alpha\alpha}$ is invertible for any $k$ and any $\alpha$. Now we have the direction. If the full Newton step size is adopted, i.e., $\mathbf{z}^{k+1}=\mathbf{z}^{k}+\mathbf{d}^{k}$, then (3.2) implies

[TABLE]

Therefore, the updated point $\mathbf{z}^{k+1}$ is feasible to the problem (1.6). However, the full Newton step size generally does not guarantee the descent property of the objective function, that is, $f(\mathbf{z}^{k+1})\leq f(\mathbf{z}^{k})$ can not be ensured. To overcome such a drawback, we exploit the following operator

[TABLE]

For some carefully chosen $\sigma_{k}$, we set $\mathbf{z}^{k+1}:=\mathbf{z}^{k}(\sigma_{k})$. Then we will show that in this way, $\mathbf{z}^{k+1}$ is not only always feasible to the problem (1.6) but also satisfies the descent property (see Lemma 3.3). Now we summarize the whole framework of Newton method in Algorithm 1.

Remark 3.1

Regarding NSLR, we have some comments.

• i)

To pick an $\alpha\in\Sigma_{s}(\mathbf{z}^{k}-\tau{\boldsymbol{g}}^{k})$, only $s$ largest elements (in absolute) are selected, which enables us to use a MATLAB built-in function mink. The computational complexity is about $\mathcal{O}(p+s\log s)$.

• ii)

Updating $\mathbf{d}^{k}$ by (3.2) involves two main calculations $H^{k}_{\alpha\alpha}$ with $|\alpha|=s$ and its inverse. Their computational complexities are $\mathcal{O}(sn+s^{2}n)$ and $\mathcal{O}(s^{3})$. So, the whole complexity is $\mathcal{O}(s^{3}+s^{2}n)$, which means the computation is quite fast if $\max\{s,n\}\ll p$.

• iii)

For a halting condition, we will calculate the quantity $\|\theta_{\alpha}^{k}\|$. If $\|\theta_{\alpha}^{k}\|=0$, then $\mathbf{z}^{k}$ satisfies (2.17). This means $\mathbf{z}^{k}$ is a $\tau$-stationary point by Theorem 2.4. Therefore, it makes sense to terminate NSLR if the quantity $\|\theta_{\alpha}^{k}\|$ is sufficiently small.

3.2 Global and quadratic convergence

To derive the convergence properties, we denote some notation hereafter. Let $\beta$ be the index set related the previous iteration $\mathbf{z}^{k-1}$ selected by

[TABLE]

Based on $\mathbf{z}^{k}=\mathbf{z}^{k-1}(\sigma_{k-1}),\mathbf{z}^{k+1}=\mathbf{z}^{k}(\sigma_{k})$ in Algorithm 1 and the definition of $\mathbf{z}^{k}(\sigma)$ in (3.3), we must have

[TABLE]

Let ${\gamma}:=\beta\setminus{\alpha}$. Then one can observe that

[TABLE]

This gives rise to the following properties

[TABLE]

Based on these, we have the following results.

Lemma 3.2

Let $\{\mathbf{z}^{k}\}$ be the sequence generated by Algorithm 1. We have the following properties:

[TABLE]

Proof It follows from (3.2) that

[TABLE]

Direct calculation yields the following chain of equations,

[TABLE]

which leads to the truth

[TABLE]

where the inequality is from the fact that $f$ being strongly convex with the constant $\lambda$ and strongly smooth with the constant $L$ so as to satisfy

[TABLE]

For the part $\overline{\alpha}$, since $\alpha\in\Sigma_{s}(\mathbf{z}^{k}-\tau{\boldsymbol{g}}^{k})$, the definition of $\alpha\in\Sigma_{s}$ in (2.16) implies

[TABLE]

Now for $i\in\alpha\setminus\beta=:\eta$, we have $z_{i}^{k}=0$ due to (3.4). Then the above condition and ${\gamma}\subseteq\overline{\alpha}$ result in

[TABLE]

which together with $\mathbf{z}_{{\overline{\alpha}}}^{k}=-\mathbf{d}_{{\overline{\alpha}}}^{k}$ from (3.2) and $|\eta|=|\alpha|-|\alpha\cap\beta|=s-|\alpha\cap\beta|=|{\beta}|-|\alpha\cap\beta|=|{\gamma}|$ leads to

[TABLE]

The above condition allows us to derive that

[TABLE]

where the last inequality is owing to (3.12) and

[TABLE]

showing the desired results. \qed

The first result below shows that the Newton direction $\mathbf{d}^{k}$ is a descent direction and the Amijio-type step size $\sigma_{k}$ always exists and is away from zero.

Lemma 3.3** (Descent property)**

Let $\{\mathbf{z}^{k}\}$ be the sequence generated by Algorithm 1 and

[TABLE]

Denote $C:=\min\left\{{1}/{\tau}-L-\tau L^{2},~{}2\lambda-\tau L^{2}\right\}>0$. Then we have

[TABLE]

Moreover, for any $c\overline{\sigma}\leq\sigma\leq\overline{\sigma}:=\lambda/(2L)$, it holds

[TABLE]

This indicates $\inf_{k\geq 0}\sigma_{k}\geq c\overline{\sigma}>0$.

Proof  Denote two parameters

[TABLE]

It is easy to see that $C_{2}>0$ by (3.24) and $C_{1}>0$ due to

[TABLE]

and $0<c,\lambda/L<1$. Overall, $C>0$. Then by (3.15), we have

[TABLE]

Note from (3.3) that

[TABLE]

The strong smoothness of $f$ with the constant $L$ yields

[TABLE]

We next show $F\leq 0$. It follows from (3.15) that

[TABLE]

One can check that, by $0<\tau<\overline{\tau}\leq{c\lambda^{2}}/({4L^{3}})$, it follows

[TABLE]

where the last inequality is true if $c\lambda/(2L)\leq\sigma\leq\lambda/(2L)$. For the $\overline{\alpha}$ part,

[TABLE]

where the above three inequalities used the facts that

• (a)

$\tau L^{2}\leq{c\lambda^{2}}/({4L})\leq\lambda(c/4)(\lambda/{L})\leq L$;

• (b)

$\sigma\leq 1$ due to $\sigma_{k}=c^{r}$ and $c\in(0,1)$ in Algorithm 1;

• (c)

$\tau\leq{c\lambda^{2}}/({4L^{3}})\leq 1/({4L})$.

Therefore, $F\leq 0$, displaying (3.26). Then the Armijo step-size rule indicates that $\inf_{k\geq 0}\sigma_{k}\geq{c\lambda}/({2L})>0$. \qed

Now we are ready to display the main results of the method including the global convergence to a $\tau$-stationary point, the support set identification, and the quadratic convergence rate.

Theorem 3.4** (Global and quadratic convergence)**

Let $\{\mathbf{z}^{k}\}$ be the sequence generated by Algorithm 1 and $0<\tau<\overline{\tau}$. The following results hold.

• i)

The whole sequence converges to a $\tau$-stationary point denoted by $\mathbf{z}^{*}$, which is at least a local minimizer of the problem (1.6).

• ii)

For sufficiently large $k$, the support sets of the sequence are identified by

[TABLE]

• iii)

The sequence converges to $\mathbf{z}^{*}$ quadratically, namely,

[TABLE]

Proof i) Lemma 3.3 shows that $\sigma_{k}\geq c\overline{\sigma}$ and results in

[TABLE]

Then it follows from the above inequality that

[TABLE]

where the last inequality is due to $f$ is positive. Hence $\|\mathbf{d}^{k}\|\rightarrow 0$, which suffices to $\|\mathbf{z}^{k+1}-\mathbf{z}^{k}\|\rightarrow 0$ since

[TABLE]

This also indicates $\|\mathbf{z}^{k}_{\overline{\alpha}}\|^{2}=\|\mathbf{d}^{k}_{\overline{\alpha}}\|^{2}\rightarrow 0$ and by (3.23) suffices to

[TABLE]

Let $\{\mathbf{z}^{k_{\ell}}\}$ be the convergent subsequence of $\{\mathbf{z}^{k}\}$ that converges to $\mathbf{z}^{*}$ and

[TABLE]

Since there are only finitely many choices for $\alpha_{\ell}\subseteq[p]$, (re-subsequencing if necessary) we may without loss of any generality assume that the sequence of the index sets $\{\alpha_{\ell}\}$ shares a same index set, denoted as $\alpha_{\infty}$. That is

[TABLE]

Now by letting ${\boldsymbol{g}}^{*}:=\nabla f(\mathbf{z}^{*})$, one can show that

[TABLE]

In addition, the definition of $\Sigma_{s}$ in (2.16) implies

[TABLE]

Taking the limit of both sides of the above inequality yields

[TABLE]

Here, we used the facts that (3.39) and $\|\mathbf{z}^{k_{\ell}}_{\overline{\alpha}_{\infty}}\|=\|\mathbf{z}^{k_{\ell}}_{\overline{\alpha}_{\ell}}\|\rightarrow 0$. Since $\mathbf{z}^{k_{\ell}}\rightarrow\mathbf{z}^{*}$, we have $\mathrm{supp}(\mathbf{z}^{*})\subseteq\alpha_{\infty}$.

If $\|\mathbf{z}^{*}\|_{0}=s$, then $\mathrm{supp}(\mathbf{z}^{*})=\alpha_{\infty}$, which by (3.40) derives $|\tau g^{*}_{j}|\leq[\mathbf{z}^{*}]^{\downarrow}_{s},j\notin\mathrm{supp}(\mathbf{z}^{*})$. This together with (3.39) results in condition (2.8).

• 2.

If $\|\mathbf{z}^{*}\|_{0}<s$, then $\mathrm{supp}(\mathbf{z}^{*})\subset\alpha_{\infty}$, which by (3.40) delivers $|\tau g^{*}_{j}|\leq[\mathbf{z}^{*}]^{\downarrow}_{s}=0,j\notin\overline{\alpha}_{\infty}$. This by (3.39) yields ${\boldsymbol{g}}^{*}=0$, which also satisfies (2.8).

Overall, both cases show $\mathbf{z}^{*}$ is a $\tau$ stationary point of (1.6). From Theorem 2.3, a $\tau$ stationary point of (1.6) is a unique local minimizer, namely, $\mathbf{z}^{*}$ is isolated. By [41, Lemma 4.10], the whole sequence converges to $\mathbf{z}^{*}$ because $\mathbf{z}^{*}$ is isolated and $\|\mathbf{z}^{k+1}-\mathbf{z}^{k}\|\rightarrow 0$.

ii) We proved that the whole sequence converges to $\mathbf{z}^{*}$. Denote $\alpha_{*}:=\mathrm{supp}(\mathbf{z}^{*})$. If $\mathbf{z}^{*}=0$, then the conclusion holds clearly due to $\alpha_{*}=\emptyset$. Consider $\mathbf{z}^{*}\neq 0$. For sufficiently large $k$, we must have

[TABLE]

If $\alpha_{*}\nsubseteq\mathrm{supp}(\mathbf{z}^{k})$, then there is an $i_{0}\in\alpha_{*}\setminus\mathrm{supp}(\mathbf{z}^{k})$ satisfying

[TABLE]

which is a contradiction. Therefore, $\alpha_{*}\subseteq\mathrm{supp}(\mathbf{z}^{k})$. By (3.4), we have $\mathrm{supp}(\mathbf{z}^{k})\subseteq\beta$, where $|\beta|=s$ by (2.16). Therefore, if $\|\mathbf{z}^{*}\|_{0}=s$ then $\beta\equiv\mathrm{supp}(\mathbf{z}^{k})\equiv\alpha_{*}$ for any sufficiently large $k$. Particularly, $\alpha_{*}\equiv\mathrm{supp}(\mathbf{z}^{k+1})\equiv\alpha$. If $\|\mathbf{z}^{*}\|_{0}<s$ then $\alpha_{*}\subseteq\mathrm{supp}(\mathbf{z}^{k})$ and $\alpha_{*}\subseteq\mathrm{supp}(\mathbf{z}^{k+1})\subseteq\alpha$ from (3.4). So (3.35) is true.

iii) For sufficiently large $k$, it follows from $\alpha_{*}\subseteq\alpha$ by ii) that $\mathbf{z}^{*}_{\overline{\alpha}}=0$. Since $\mathbf{z}^{*}$ is a $\tau$-stationary point, (2.8) indicates ${\boldsymbol{g}}^{*}=0$ if $\|\mathbf{z}^{*}\|_{0}<s$ and ${\boldsymbol{g}}^{*}_{\alpha_{*}}=0$ if $\|\mathbf{z}^{*}\|_{0}=s$. While for the latter case, there is $\alpha_{*}=\alpha$ by ii). Overall, we have

[TABLE]

For any $0\leq t\leq 1$, denote $\mathbf{z}(t):=\mathbf{z}^{*}+t(\mathbf{z}^{k}-\mathbf{z}^{*})$ and $H^{k}(t):=\nabla^{2}f(\mathbf{z}^{k}(t))$. Then (2.4) derives

[TABLE]

Moreover, by Taylor expansion, one has

[TABLE]

From (3.4) and (3.41), we have $\mathbf{z}^{k+1}_{\overline{\alpha}}=\mathbf{z}^{*}_{\overline{\alpha}}=0$ and the following relations

[TABLE]

where the first inequality is due to $\|\cdot\|^{2}$ being convex and the last inequality is from Lemma 3.3 that $c\overline{\sigma}\leq\sigma_{k}\leq\overline{\sigma}$. For the second term in (3.48), we have

[TABLE]

where the forth inequality used a fact that

[TABLE]

It follows from $\mathbf{d}^{k}_{\overline{\alpha}}=-\mathbf{z}^{k}_{\overline{\alpha}}$ and (3.41) that

[TABLE]

leading to the following fact

[TABLE]

Now we have three facts: (3.60), $\mathbf{z}^{k}\rightarrow\mathbf{z}^{*}$ from i), and $\langle{\boldsymbol{g}}^{k},\mathbf{d}^{k}\rangle\leq-C\|\mathbf{d}^{k}\|^{2}$ from Lemma 3.3. They and [42, Theorem 3.3] allow us to claim that eventually the step size $\sigma_{k}$ determined by the Armijo rule is 1, namely, $\sigma_{k}=1$. Then it follows from (3.48) that

[TABLE]

delivering the quadratic convergence property of the sequence. \qed

4 Numerical Experiments

In this part, we will conduct extensive numerical experiments of NSLR111Available at https://github.com/ShenglongZhou/NSLR by using MATLAB (R2017b) on a desktop of 8GB of memory and Inter Core i5 2.7Ghz CPU, against seven leading solvers on both synthetic and real datasets.

4.1 Test examples

We first adopt two types of randomly generated data: the one with the i.i.d. features $(\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n})\in\mathbb{R}^{p\times n}$ and the one with independent features with each of $\mathbf{x}_{i}$ being generated by an autoregressive process [43]. Then eight real datasets are taken into consideration to test the selected methods.

Example 4.1** (Independent Data [36, 38])**

To generate data labels $\mathbf{y}\in\{0,1\}^{n}$, we first randomly divide $[n]$ into two parts and set $y_{i}=0$ for one part and $y_{i}=1$ for the other. Then the feature data is produced by

[TABLE]

with $\mathbb{R}\ni v_{i}\sim\mathcal{N}(0,1)$, $\mathbb{R}^{p}\ni\mathbf{w}_{i}\sim\mathcal{N}(0,I_{p})$, where $\mathcal{N}(0,{I}_{p})$ is the normal distribution with mean zero and variance identity. Here, 1 is the vector with all entries being ones. Since the sparse parameter $\mathbf{z}^{*}\in\mathbb{R}^{p}$ is unknown, different $s(<n)$ will be tested to pursue a sparse solution.

Example 4.2** (Correlated Data [44, 2])**

The sparse parameter $\mathbf{z}^{*}\in\mathbb{R}^{p}$ has $s$ nonzero entries drawn independently from the standard Gaussian distribution. Each data sample $\mathbf{x}_{i}=(x_{i1},x_{i2},\ldots,x_{ip})^{T},i\in[n]$ is an independent instance of the random vector generated by an autoregressive process [43] determined by

[TABLE]

with $x_{i1}\sim\mathcal{N}(0,1)$, $v_{ij}\sim\mathcal{N}(0,1)$ and $\rho\in[0,1]$ being the correlation parameter. The data labels $y_{i}\in\{0,1\}$ are then drawn randomly according to the Bernoulli distribution with the conditional probability (1.1).

Example 4.3** (Real data)**

Eight real data sets are taken into consideration. They are summarized in Table 1, where arcene and newsgroup are taken from UCI repository222http://archive.ics.uci.edu/ml/index.php and glmnet package333https://web.stanford.edu/~hastie/glmnet_matlab/, and the rest of them are LIBSVM data444https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/. Moreover, all datasets are feature-wisely scaled to $[-1,1]$. All $-1$s in the label classes $\mathbf{y}$ are replaced by 0. The sizes of training data and testing data are denoted by $m_{1}$ and $m_{2}$ respectively.

4.2 Implementation

In the model (1.6), we set $\lambda=10^{-5}/n$. As mentioned in Remark 3.1, we terminate the method if $\|\theta^{k}\|<10^{-10}\sqrt{p}$ or $k>2000$. For the starting point and parameters of NSLR, we set $\mathbf{z}^{0}=0$ and $c=0.5$. For the parameter $\tau$, Theorem 3.4 indicates $0<\tau<\overline{\tau}$. While this is a sufficient condition. So it is unnecessary to choose a $\tau$ to satisfy the condition strictly, not to mention, $\overline{\tau}$ might be too small if we set a tiny $\lambda$.

Alternative is to pick a proper $\tau$ by tuning it from a wide range of values. For instance, we tested NSLR on a range of selections of $\tau=10^{\varrho}$ with $\varrho\in\{-3,-2.8,\ldots,1,1.2,1.4\}$ for solving Example 4.1. As reported in Fig. 1, for $s=5$, NSLR generated the best results when $\tau$ was around $3.98$, while for $s=25$, it delivered the best results when $\tau$ was around $10$. Therefore, for different scenarios, the best option $\tau$ may be varied, which indicates the manual selection of $\tau$ is necessary to achieve a better performance. This apparently would incur expensive computational costs.

However, the empirical numerical experience have demonstrated that adaptively updating parameters during the process is a practicable strategy. Hence, we start $\tau$ with a fixed scalar $\tau_{0}=15$ and update $\tau_{k+1}=0.75\tau_{k}$ if $k$ is a multiple of 10 and $\|\theta^{k}\|>1/k$, and $\tau_{k+1}=\tau_{k}$ otherwise. In the sequel, we adopt this strategy for NSLR and the numerical comparisons with other leading solvers will show the superior performance of our method under such a strategy.

4.3 Benchmark methods

Since there is an impressive body of methods that have been developed to address the sparse logistic regression, we only focus on those programmed by Matlab. Solvers with codes being online unavailable or being written by other languages, such as R and C, are not selected for comparisons. We thus choose 7 solvers mentioned in Section 1.2, which should be enough to make comprehensive comparisons. We summarize them into the following table.

For SLEP, we use it to solve (1.4) with $\phi_{\nu}(\mathbf{z})=\nu_{1}\|\mathbf{z}\|^{2}_{2}+\nu_{2}\|\mathbf{z}\|_{1}$, whilst IRLS-LARS aims to solve the case of $\nu=0$. APG and GIST are taken to solve the capped $\ell_{1}$ logistic regression with $\phi_{\nu}(\mathbf{z})=\nu_{3}\min(|\mathbf{z}_{i}|,\nu_{4})$. We only use non-monotonous version of APG since its numerical performance was better than that of the monotonous version [30]. Note that methods that aim at solving the model (1.4) involve a penalty parameter $\nu$, whilst those tackling (1.6) need the sparsity level $s$. To make results comparable, we adjust their default parameters $\nu$ for each method to guarantee the generated solution $\mathbf{z}$ satisfying $\|\mathbf{z}\|_{0}\leq p/2$. We will report three indicators: $(\ell(\mathbf{z}),~{}\texttt{SER},~{}\texttt{Time})$ to illustrate the performance of methods, where Time (in seconds) is the CPU time, $\mathbf{z}$ is the solution obtained by each method and SER is the sign error rate defined by

[TABLE]

Here $\mathrm{sign}(a_{+})$ is the sign of the projection of $a$ onto a non-negative space, namely, it returns $1$ if $a>0$ and [math] otherwise.

4.4 Numerical comparison

We now report the performance of eight methods on the above three examples. To avoid randomness, we report average results over 10-time independent trails for Examples 4.1 and 4.2 since they involve in randomly generated data.

(a) Comparison on Example 4.1. To observe the influence of the sparsity level $s$ on four greedy methods: NSLR, GPGN, GraSP and NTGP, we fix $p=10000,n=p/5$ and alter $s\in\{400,600,\ldots,1600\}.$ As demonstrated in Fig. 2, NSLR outperforms others in terms of the lowest $\ell(\mathbf{z})$ and SER and the shortest time, followed by GPGN. By contrast, GraSP always performs the worst results, which means this first-order method is not competitive when against the other three methods, three second-order methods.

To observe the influence of the ratio of the sample size $n$ and the number of features $p$ on all eight methods, we fix $p=20000,\ s=0.1p$ and vary $n/p\in\{0.1,0.2,\ldots,0.7\}$. Apart from recording the four indicators, we also report the number of non-zeros of the solution $\mathbf{z}$ generated by each method. Here, for the LARS, we stop it when $s$ variables are selected and their default stopping conditions are met since LARS only adds one variable at each iteration (see [45]). We set $\nu_{1}=10^{-1},\nu_{2}=10^{-2}$ for SLEP, $\nu_{3}=10^{-2},\ \nu_{4}=10^{-4}$ for APG and $\nu_{3}=10^{-3}{\tt abs(randn)},\ \nu_{4}=10^{-5}{\tt abs(randn)}$ for GIST.

As presented in Fig. 3, in terms of $\ell(\mathbf{z})$ and SER, again NSLR performs the best results, followed by GPGN and GIST. It is obvious that LARS and SLEP produce undesirable results compared with other methods. For the computational time, NSLR runs the fastest, while GraSP  and APG run relatively slow with over 1000 seconds when $n/p\geq 0.6$. Table 3 shows the sparsity levels $\|\mathbf{z}\|_{0}$ only in LARS is lower than our NSLR. This is because LARS fails to recover the support and vanishes when $s=500$ in this numerical experiment (this phenomenon had also been observed in [45, 46].)

(b) Comparison on Example 4.2. To observe the influence of the correlation parameter $\rho$ on eight methods, we set $p=1000,s=0.1p$ but choose $\rho\in\{0,1/3,1/2\}$. Fig. 4 shows the average $\ell(\mathbf{z})$ gotten by eight methods for a wide range of the ratio $n/p\in\{0.1,0.2,\ldots,0.9\}$. Apparently, at three different values of $\rho$, NSLR always performs stably best results. Moreover, the trends in these eight methods perform generally consistent, which indicates the correlation parameter has little influence on these methods. Therefore, we further fix $\rho=1/2$ and observe the performance of eight methods under higher dimensions.

Now we alter $p\in\{10000,20000,30000\}$ with $n=0.2p,s=0.05p$ or $s=0.1p$. Here, we set $\nu_{1}=10^{-2},\nu_{2}=10^{-1}$ for SLEP  $\nu_{3}=10^{-2},\ \nu_{4}=5\times 10^{-4}$ for APG and $\nu_{3}=5\times 10^{-3}{\tt abs(randn)},\ \nu_{4}=5\times 10^{-5}{\tt abs(randn)}$ for GIST. As reported in Table 4, for cases of $p=20000,30000$, LARS basically fails to render desirable solutions due to the highest $\ell(\mathbf{z})$. It can be clearly seen that NSLR always provides the best accuracies with consuming shortest time.

(c) Comparison on Example 4.3. To observe the performance for above all eight methods on real data sets, we select eight real data sets with different dimensions. The highest dimension is up to millions (see news20.binary). Table 5 reports results for eight methods on four datasets without testing data ($m_{2}=0$): arcene, colon-cancer, news20.binary and newsgroup. For the last two datasets, LARS makes our desktop run out of memory, thus its results are omitted here. Clearly, NSLR is more efficient than others for all test instances. For example, NSLR only uses 1.365 seconds for data newsgroup with $p=777811$ features and achieves the smallest logistic loss with the sparsest solution.

When all methods solve the datasets with testing data ($m_{2}>0$): duke breast-cancer, leukemia, gisette and rcv1.binary, the table is a little different. Since the testing data is taken into consideration, we add two indicators to illustrate the performance of each method: $\ell(\mathbf{z})$-test and SER-test. Results are reported in Table 6, where $\ell(\mathbf{z})$ and $\ell(\mathbf{z})$-test denote the objective function value on training data and testing data respectively, and similar to SER-train and SER-test. For cases of duke breast-cancer and leukemia, LARS stops when the maximum number of iterations reaches the $\min\{m_{1},p\}$. Hence its produced $||\mathbf{z}||_{0}$s are less than other methods. We can see that NSLR could guarantee a good performance on the testing data as well as the training data.

5 Conclusion

Despite the NP-hardness of the sparsity constrained logistic regression (1.6), we benefited from its nice properties of the objective function and introduced the $\tau$-stationary point as an optimality condition. This can be converted to an equation system that makes the Newton method effective. Since an order $s$ (which is far smaller than $p$) principal sub-matrix of the whole Hessian of the objective function is taken into account in each step, the proposed method NSLR has a relatively low computational complexity. The success in acquiring the global convergence stems from the realization of the Armijo-type line search in the method. What is more, the generated sequence also converges to a $\tau$-stationary point quadratically, which well justifies the outstanding performance of NSLR theoretically. It is worth mentioning that we reasonably extended the classical Newton method for solving unconstrained and continuous problems to the sparsity constrained logistic regression. The numerical performance against several state-of-the-art methods demonstrated that NSLR is remarkably efficient and competitive, especially in large scale settings.

We feel that the proposed method might be capable of solving the general strong convex optimization problems with the sparsity constraint. This deserves exploring in future.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (11971052) and Beijing Natural Science Foundation (Z190002). We particularly thank the referee and the Principal Editor who offered us valuable suggestions to improve this paper greatly.