上海理工大学理学院，上海

收稿日期：2022年9月18日；录用日期：2022年10月8日；发布日期：2022年10月18日

摘要

Minimax Concave Penalty (MCP)正则优化问题在诸多科学领域有着广泛的应用，例如：机器学习、信号、图像恢复以及逻辑回归等问题。本文基于MCP正则项研究了一种二阶加速优化算法，该算法的主要思想是在于如何得到一个使得目标快速下降的方向，主要的方法是设计对偶半光滑牛顿法求解子问题。我们基于所提出的模型提出新的算法。通过数值试验部分的对比验证了我们所提出的算法的有效性和高效性。

关键词

非凸优化，临近算子，半光滑牛顿法

Efficient Second-Order Algorithm for MCP Regular Optimization Problems

Yuanhao Zhang

College of Science, University of Shanghai for Science and Technology, Shanghai

Received: Sep. 18^th, 2022; accepted: Oct. 8^th, 2022; published: Oct. 18^th, 2022

ABSTRACT

The Minimax Concave Penalty (MCP) regularization optimization problem is widely used in many scientific fields, such as machine learning, signal restoration, image restoration and logistic regression. This paper studies a second-order accelerated optimization algorithm based on MCP regularization term. The main idea of the algorithm is how to get a direction that makes the target drop quickly. The main method is to design dual semi-smooth Newton method to solve the subproblem. We propose a new algorithm based on the proposed model. The effectiveness and efficiency of the proposed algorithm are verified by numerical experiments.

Keywords:Non-Convex Optimization, Proximal Operator, Semi-Smooth Newton Method

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

本文考虑一类如下的优化问题：

${\min }_{x\in R^n}F\left(x\right):=f\left(x\right)+\eta \hat{g}\left(x\right)$ (1-1)

其中f是光滑的并且梯度是Lipschitz连续的， $\hat{g}$ 是正常的下半连续函数(可能非凸非光滑)， $\eta >0$ 是给定参数。对于加式类型的复合优化问题，近些年涌现出了许多优秀的优化方法如：临近梯度方法(PGA) [1]、交替方向乘子法(ADMM) [2]、DC临近牛顿法 [3]、半光滑牛顿法 [4] 等。并且该类型的优化问题在实际生活中有着广泛的应用，例如：非线性最小二乘 [5]、鲁棒相位恢复 [6]、低秩矩阵补全 [7]、信号恢复 [8]、逻辑回归 [9] 以及诸多机器学习问题。而随着科学技术的进一步发展，以及大数据时代的到来，优化算法在实际生活中的应用也是越来越广泛，从简单的凸的无约束优化算法到现在迅速发展的非凸带约束的优化算法，表明了人类永不止步的探索精神。从简单到难也意味着我们对于数据的处理以及要求越来越高，对应用于实际问题的算法要求越来越严格。随着时间的推移对于凸优化问题的研究已经发展的相当成熟，而对于非凸优化问题的研究还是寥寥无几，所以近几年来有大批学者着手于发展研究带有非凸非光滑正则项的优化问题，而大多数的非凸优化都是以求解稀疏解为目的的。就稀疏优化问题而言，最初的探究想法是将函数 $\hat{g}$ 设置为 $L_0$ 范数，但是由于带 $L_0$ 范数的优化模型是一个NP难的问题，而目前大部分的解决方法是利用 $L_1$ 范数对其进行一个凸的替代，这样非凸的问题就变成了一个凸的优化问题，我们就可以使用现有的各种高效算法来求解，但事实上在实际生活中的许多工业或者自然科学等研究问题中目标大多是一个非凸的函数，这样我们原本的凸优化方法就显得不那么高效，由此非凸优化算法也就应运而生了。非凸优化问题在上述应用中也都有所发展。本文研究带非凸正则项的优化问题，目前研究发现一些性质比较好的非凸惩罚有Minimax Concave Penalty (MCP) [2] [8] [10] [11] [12] [13] [14]、Smoothly Clipped Absolute Deviation (SCAD) [12] [15] 等。最近的一些研究中发现了基于上述非凸正则项的一些变种同样有着良好的发展潜力有待我们研究，例如：GMCP [16] [17]，VMCP [13] 等等。尽管对于现阶段的各类非凸惩罚已经有了大量的研究，但是大多数都是一阶方法，例如：邻近梯度法(PGA) [8] [18]、交替方向乘子法(ADMM) [2] [15]、坐标下降法(CD) [19] 等。对于带非凸正则项的优化问题，二阶算法的研究就显得相对较少。

随着大数据时代的发展，显然凸优化方法一定程度上已经无法满足自然科学的更高要求，经过一段时间的研究和发现，非凸优化问题的求解虽然困难，但是通过变换也可以有效地解决，重要的是非凸优化问题能够弥补凸优化问题中的一些不足，非凸优化问题相比于凸优化问题而言优秀的算法虽然较少，但是目前仍有大批国内外学者在竭尽全力研究。其中较为突出的是非凸稀疏正则优化问题，其中性质良好的MCP正则项2010年由Zhang [14] 在理论和算法上作出了巨大的贡献，随后2015年由Li和Lin [18] 提出了解决非凸非光滑优化问题的加速近端梯度算法(APG)，其中单调APG具体框架如下：

$\begin{array}{l}{y}_k=x_k+\frac{t_{k-1}}{t_k}\left(z_k-x_k\right)+\frac{t_{k-1}-1}{t_k}\left(x_k-x_{k-1}\right)\\ z_{k+1}=pro{x}_{{\alpha }_{Y}g}\left(y_k-{\alpha }_Y\nabla f\left(y_k\right)\right)\\ v_{k+1}=pro{x}_{{\alpha }_{X}g}\left(x_k-{\alpha }_X\nabla f\left(x_k\right)\right)\\ t_{k+1}=\frac{\sqrt{4{\left(t_k\right)}^2+1}+1}{2}\\ x_{k+1}=\{\begin{array}{l}{z}_{k+1},\text{if}F\left(z_{k+1}\right)\le F\left(v_{k+1}\right)\\ v_{k+1},\text{otherwise}\end{array}\end{array}$

之后，2019年Shi和Huang [12] 提出了一种二阶的，针对高维非凸稀疏优化问题的半光滑牛顿算法，并证明了他们所提出的算法局部超线性收敛于KKT点，同时数值试验也验证了该方法计算效率更高。但是，尽管该方法是二阶方法收敛速度快，而文中仅仅只考虑了残差下降并没有考虑目标的下降，其次该文章也只证明了局部的超线性收敛速度，并没有考虑全局收敛性。

在2021年Wu和Li [8] 等人推广了近端梯度算法，把传统近端项中的欧式范数改为了更一般的Bregman范数，进而设计出了新的一阶算法，所提出的新算法一般框架如下：

$\begin{array}{l}{y}^k:=x^k+{\alpha }^k\left(x^k-x^{k-1}\right)\\ z^k:=x^k+{\beta }^k\left(x^k-x^{k-1}\right)\\ x^{k+1}\in {\mathcal{P}}_{{\lambda }_{k}g}^h\left(x^k;\nabla f\left(z^k\right)-\frac{1}{{\lambda }_k}\left(y^k-x^k\right)\right)\end{array}$

然后在KL性质的假设条件下证明了全局收敛性和收敛速度，同时以SCAD和MCP两种非凸正则项为例通过数值实验，验证了该算法的有效性和高效性。但是由于他们所提出的算法仍属于一阶算法，收敛速度依然有待提高。

2. 相关基础知识

基础知识

为了方便后续文章模型的建立以及算法的构建，首先我们明确一些假设、符号、定义以及方法(对于基本的符号、定义和方法可以参考 [11] [20] [21])。

定义2.1. (梯度)给定函数 $f:R^n\to R$，且f在点x的一个邻域内有意义，若存在向量 $\stackrel{\dot{}}{g}\in R^n$ 满足：

$\underset{\delta \to 0}{\lim }\frac{f\left(x+\delta \right)-f\left(x\right)-{\stackrel{\dot{}}{g}}^{\text{T}}\delta }{\Vert \delta \Vert }=0$

就称f在点x处可微。此时 $\stackrel{\dot{}}{g}$ 称为f在点x处的梯度。并记任意一阶可微函数 $f\left(x\right)$ 梯度为： $\nabla f\left(x\right)$。

假设2.2. 假设目标函数 $F\left(x\right)$ 有下界，即 $F\left(x\right)\ge \bar{F}$，函数 $f:R^n\to R$ 是二次连续可微的，其梯度是Lipschitz连续的。对 $\forall k$ 满足 $\left[{\lambda }_{\min }\left(B_k\right)-\eta \delta \right]\underset{\_}{\succ }0$。 $\hat{g}$ 为MCP正则项。

其中函数f的梯度为Lipschitz常数为： $L_{\nabla f}$，$B={\nabla }^{2}f$，$B_k$ 代表第k次迭代对 ${\nabla }^{2}f$ 的近似。 ${\lambda }_{\min }\left(A\right)$ 是矩阵A的特征值的最小值。 $\eta$ 是罚参数， $\delta$ 是与正则项相关的参数。

对于任意实值扩张函数 $f:R^n\to \left(-\infty ,+\infty \right]$，它的定义域如下：

$\text{dom}f:=\left\{x\in R^n:f\left(x\right)<+\infty \right\}$ (2-1)

接着我们给出临近映射和Moreau-envelope的定义分别如下：

$pro{x}_{\mu f}\left(x\right):=\arg {\min }_u\left\{f\left(u\right)+\frac{1}{2\mu }{\Vert u-x\Vert }^2\right\}$ (2-2)

$M_f^{\mu }\left(x\right):={\inf }_u\left\{f\left(u\right)+\frac{1}{2\mu }{\Vert u-x\Vert }^2\right\}$ (2-3)

欧式空间上的内积和范数定义如下：

$\begin{array}{l}\langle x,y\rangle =x^{\text{T}}y\\ {\Vert x\Vert }_2=\sqrt{\langle x,x\rangle }\end{array}$

对于正则项 $\hat{g}\left(x\right)$，我们根据文献 [14] 的研究，采用性质较好和研究较多的非凸非光滑正则项MCP，其具体定义如下：

$\hat{g}\left(x\right)={\mathcal{P}}_{MC}\left(x;\lambda ,\gamma \right)={\displaystyle {\sum }_{i=1}^n\left\{P_{MC}\left(x_i;\lambda ,\gamma \right)\right\}}$ (2-4)

其中 $x\in R^n$，$P_{MC}\left(x_i;\lambda ,\gamma \right)$ 定义如下：

$P_{MC}\left(x_i;\lambda ,\gamma \right)=\lambda {\displaystyle {\int }_0^{\left|x_i\right|}{\left(1-z/\left(\gamma \lambda \right)\right)}_{+}\text{d}z}=\{\begin{array}{l}\lambda \left|x_i\right|-\frac{1}{2\gamma }{x}_i^2,\text{if}\left|x_i\right|\le \lambda \gamma \\ \frac{1}{2}\gamma {\lambda }^2,\text{if}\left|x_i\right|>\lambda \gamma \end{array}$ (2-5)

根据定义我们容易得到函数 $\alpha \hat{g}\left(x\right)$ (其中 $\alpha \in \left(0,\gamma \right)$ )的临近映射为： $\text{firm}\left(x;\alpha \lambda ,\lambda \gamma \right)$ [13]，具体表达式如下：

$\begin{array}{c}pro{x}_{\alpha \hat{g}}\left(x\right)=\arg {\min }_u\left\{\hat{g}\left(u\right)+\frac{1}{2\alpha }{\Vert x-u\Vert }^2\right\}\\ =\text{firm}\left(x;\alpha \lambda ,\lambda \gamma \right)\\ =\{\begin{array}{ll}0,\hfill & \text{if}\left|x\right|\le \alpha \lambda \hfill \\ \text{sign}\left(x\right)\frac{\gamma \left(\left|x\right|-\alpha \lambda \right)}{\gamma -\alpha },\hfill & \text{if}\alpha \lambda <\left|x\right|\le \lambda \gamma \hfill \\ x,\hfill & \text{if}\left|x\right|>\lambda \gamma \hfill \end{array}\end{array}$

进而根据文献 [4] 知道函数 $pro{x}_{\alpha \hat{g}}\left(x\right)$ 是广义可微的，并且容易得到 $pro{x}_{\alpha \hat{g}}\left(x\right)$ 的广义雅可比是一个对角阵，其对角元素为：

$J_{mcp}^{ii}\left(x\right)=\{\begin{array}{ll}0,\hfill & \text{if}\left|x_i\right|<\alpha \lambda \hfill \\ \left[0,\frac{\gamma }{\gamma -\alpha }\right],\hfill & \text{if}\left|x_i\right|=\alpha \lambda \hfill \\ \frac{\gamma }{\gamma -\alpha },\hfill & \text{if}\alpha \lambda <\left|x_i\right|<\lambda \gamma \hfill \\ \left[1,\frac{\gamma }{\gamma -\alpha }\right],\hfill & \text{if}\left|x_i\right|=\lambda \gamma \hfill \\ 1,\hfill & \text{if}\left|x_i\right|>\lambda \gamma \hfill \end{array}$ (2-6)

3. 模型、算法及实验结果

3.1. 模型提出

在有了基本的符号和知识后，本节我们将给出文章具体考虑的模型以及所提出的算法，最终通过数值实验来验证我们算法的有效性以及高效性。本着重于求解一个使得目标快速下降的方向，然后对于如何求解该方向而设计快速算法。首先我们给出子问题在满足假设2.2的情况下，若当前迭代点为 $x_k$，我们假设下一步迭代点为 $x_k+d$，其中d是迭代搜索方向，将函数f在 $x_k+d$ 这点进行二阶近似展开得到：

$f\left(x_k+d\right)\approx f\left(x_k\right)+\nabla f{\left(x_k\right)}^{\text{T}}d+\frac{1}{2}{d}^{\text{T}}{B}_{k}d$ (3-1)

我们令： $\phi \left(d\right):=\nabla f{\left(x_k\right)}^{\text{T}}d+0.5d^{\text{T}}{B}_{k}d+\eta \hat{g}\left(x_k+d\right)$，进一步得到子问题形式如下：

${\min }_d\phi \left(d\right):=\nabla f{\left(x_k\right)}^{\text{T}}d+\frac{1}{2}{d}^{\text{T}}{B}_{k}d+\eta \hat{g}\left(x_k+d\right)$ (3-2)

问题(3-2)就是本文主要考虑的模型。根据文献 [13] 我们知道，当 $\delta \ge \frac{1}{\gamma }$ 时， $\hat{g}\left(x\right)+\frac{\delta }{2}{\Vert y-x\Vert }_2^2$ 是一个凸函数，基于这点我们将问题(3-2)等价变形为：

$\begin{array}{l}\underset{d,y}{\min }\nabla f{\left(x_k\right)}^{\text{T}}d+\frac{1}{2}{d}^{\text{T}}\left(B_k-\eta \delta I\right)d+\eta \hat{g}\left(y\right)+\frac{\eta \delta }{2}{\Vert y-x_k\Vert }_2^2\\ \text{s}\text{.t}\text{.}{x}_k+d=y\end{array}$ (3-3)

进而我们容易得到问题(3-3)的对偶问题为：

${\max }_{\mu }{c}_{\delta }\left(\mu \right)$ (3-4)

其中 $\mu$ 为对偶问题的自变量， $c_{\delta }\left(\mu \right)$ 为：

$\begin{array}{c}{c}_{\delta }\left(\mu \right)=-\frac{1}{2}{\left[\mu -\nabla f\left(x_k\right)\right]}^{\text{T}}{\left(B_k-\eta \delta I\right)}^{-1}\left[\mu -\nabla f\left(x_k\right)\right]-\frac{{\Vert \mu \Vert }_2^2}{2\eta \delta }\\ +\eta \hat{g}\left(pro{x}_{{\delta }^{-1}\hat{g}}\left(x_k-\frac{\mu }{\eta \delta }\right)\right)+\frac{\eta \delta }{2}{\Vert x_k-\frac{\mu }{\eta \delta }-pro{x}_{{\delta }^{-1}\hat{g}}\left(x_k-\frac{\mu }{\eta \delta }\right)\Vert }_2^2\end{array}$ (3-5)

为了算法的构建我们引入一个记号： ${\Delta }_k\left(\alpha \right)\left(\alpha \in \left(0,1\right)\right)$，具体定义如下：

$\begin{array}{c}{\Delta }_k\left(\alpha \right)=\phi \left(0\right)-\phi \left(\alpha d_k\right)\\ =\eta \hat{g}\left(x_k\right)-\eta \hat{g}\left(x_k+\alpha d_k\right)-\alpha \nabla f{\left(x_k\right)}^{\text{T}}{d}_k-\frac{{\alpha }^2}{2}{d}_k^{\text{T}}{B}_k{d}_k\end{array}$ (3-6)

3.2. 算法构建

算法1. ALG1-求解问题(1-1)

算法2. 求解问题(3-4)

3.3. 数值实验及结果分析

本节我们将用随机合成的数据和实际数据来评估我们设计的算法1，为了进行比较同时我们也编码了临近梯度法(PGA)具体内容由算法3提供。从算法框架可以看出我们的算法和PGA方法的主要区别是搜索方向 $d_k$ 的计算。文章所有的代码都是用MATLAB2020a编写的，并在一台8核、芯片为Apple M1、内存为8 GB的微处理器上运行。对不同的算法将统一使用以下参数(随机数据的参数随后也会详细讨论。)：

$\epsilon ={10}^{-6},{\sigma }_1={\sigma }_2={10}^{-4},{\tau }_1={\tau }_2=0.2,\eta =0.25,\gamma =15,\lambda =0.03\ast \sqrt{2\ast \log \left(n\right)}$ (3-7)

不同的求解器停止条件我们统一使用：

$\text{NormG}:=\frac{{\Vert d_k\Vert }_2}{{\Vert d_0\Vert }_2+1}\le \epsilon$ (3-8)

算法3. PGA [20] for MCP正则优化问题

首先我们考虑随机合成的数据。实验为稀疏信号恢复问题 [8]，具体形式如下：

${\min }_{x\in R^n}\frac{1}{2}{\Vert Ax-b\Vert }_2^2+\eta \hat{g}\left(x\right)$ (3-9)

其中 $\hat{g}\left(x\right)$ 我们选择(2-4)定义的非凸正则项，矩阵 $A\in R^{m\times n}$ 是测量矩阵，向量 $b\in R^m$ 是观测信号， $\eta$ 是罚参数。我们使用MATLAB内置函数randn，sprandn，normrnd，normalize生成随机数据：

$A=\text{normalize}\left(\text{randn}\left(m,n\right),\text{'}\text{range'}\right)$

$x^{\ast }=\text{sprandn}\left(n,1,0.15\right)$

$eps=\text{normrnd}\left(0,0.01,\left[m,1\right]\right)$

$b=A{x}^{\ast }+eps$

其中 $\text{normalize}\left(\text{randn}\left(m,n\right),\text{'}\text{range'}\right)$ 是生成一个 $m\times n$ 的正态分布随机矩阵，并且按列归一化至[0, 1]。 $\text{sprandn}\left(n,1,0.15\right)$，是生成一个随机稀疏向量，非零元素约有 $n\times 1\times 0.15$ 个。 $\text{normrnd}\left(0,0.01,\left[m,1\right]\right)$ 是生成一个均值为0标准差为0.01的正态分布随机向量。我们将对 $m,n,{\sigma }_1,{\tau }_1$ 四个参数用控制变量法进行实验以验证算法的有效性。评估指标为CPU time (CPU运行时间)、Iterations (迭代次数)。图例中的eta代表参数： $\eta$，同样的gamma代表参数： $\gamma$。

首先我们设置 $m=2400$，${\sigma }_1={10}^{-4}$，${\tau }_1=0.2$，n从1000变化至1100总共11组数据。数值结果从图1中可以看出随着n增大时，相较于PGA方法，我们的算法1迭代次数有着明显的优势，在CPU运行时间方面也都远小于PGA方法。

其次我们设置 $n=1000$，${\sigma }_1={10}^{-4}$，${\tau }_1=0.2$，m从2410变化至2470总共7组数据。数值结果如图1所示，从图中可以看出我们的算法相较于PGA方法，迭代次数和CPU运行时间是远小于PGA方法的，这证明我们所提出的方法有效并且是高效的。

再者，我们将数据的维数固定在 $m=2400$，$n=1000$，$\gamma$ 从10变化至20一共6组数据，结果如图1所示，很容易看出随着 $\gamma$ 的变化我们的算法是很稳健的，无论是从迭代次数还是CPU运行时间来看，都是远远比PGA更加的优越。最后我们将以同样的维数 $m=2400$，$n=1000$，但是 ${\sigma }_1={10}^{-4}$，$\eta$ 从0.2变化至0.26一共7组数据对算法进行测试。最终的结果由图1给出。

接下来由于PGA方法运行时间较长我们将对自己的算法进行评估，评估指标使用目标函数值(Value of (f + g))和恢复精度(nx_xstar)。以不同维数和不同的 $\eta$ 值为例子说明我们算法的有效性和高效性。对于维度测试，我们将m固定在5000维，n从1000变化至3000，变化幅度为500。 $\eta =0.25$，$\gamma =15$ 其余参数不做改变。数值结果如图2所示，可以看到不同维度下我们的算法还是很稳健的。接着我们将维数保持在 $m=3000$，$n=1500$，$\eta$ 从0.22变化至0.3，变化幅度为0.02， $\gamma =15$ 其余参数不做改变。数值结果如图2所示，我们大致可以看到随 $\eta$ 值变化的时候，当 $\eta$ 保持在0.24到0.26这个区间的时候，算法效果是最优的，从而也证实了我们之前的实验设置 $\eta =0.25$ 是合适的选择。

图1. 4种参数下两算法间CPU运行时间和迭代次数的对比

图2. 算法1在不同参数下的表现

接下来我们考虑逻辑回归问题，给定数据集 $A\in R^{m\times n}$ 和标签向量 $b\in R^m$，逻辑回归问题表示如下：

$\min {\Phi }_{\eta }\left(x\right):=\frac{1}{m}{\displaystyle {\sum }_{i=1}^m\log \left(1+\exp \left(-b_i{x}^{\text{T}}{a}_i\right)\right)}+\eta \hat{g}\left(x\right)$ (3-10)

其中 $A=\left(a_1,\cdots ,a_m\right)$，$b_i\in \left\{-1,1\right\}$，$i=1,\cdots ,m$。

这部分我们将从LIBSVM [22] 数据集中选取部分实例进行测试。每个实例给出了一个矩阵A和向量b。具体的测试结果由表1和图3给出。我们将CPU运行时间(CPU(s))、算法迭代次数(Iter)和收敛精度(Res)作为本部分的评价指标，而Optval代表模型的最优值。其中收敛精度Res表示如下：

$\text{Res}=\frac{{\Vert d_k\Vert }_2}{{\Vert d_0\Vert }_2+1}\le \epsilon$

从表1可以清楚的看到，两种算法都可以有效地求解逻辑回归问题。而从图3中的两图我们可以清楚地看到，我们的算法1在20个不同的测试集上都有着优秀的表现，相比于传统的PGA方法而言迭代次数大大减少，同时也能收敛到更高的精度。虽然在表1中我们发现了测试集5我们的算法运行时间比PGA方法慢了一些，但是综合对比迭代次数和收敛的精度我我们的算法是远优于PGA方法的，所以对于时间上微小的差别基本可以忽略。从而这也验证了我们的算法在逻辑回归问题上有着卓越的表现，通过与PGA方法的对比也凸显出了我们的方法是更为有效和高效的。

综合上述实验以及结果的分析，验证了我们所设计的二阶算法在逻辑回归和稀疏信号恢复方面都是有效的，并且相较于传统的PGA方法在各个方面都更加高效，鲁棒性更好。

表1. 两算法在20组逻辑回归中的表现

图3. 20个数据集中两算法在迭代次数和收敛精度中的表现

4. 总结

本文基于带非凸正则项MCP的优化问题进行了二阶算法的研究，通过对目标可微项的二阶近似展开，得到了寻求下降方向的子问题，并利用对偶半光滑牛顿法对其进行求解。最终设计了一个新的二阶算法，通过大量的数值实验我们验证了在信号恢复和逻辑回归问题中，我们所提出的算法是有效并且高效的，而且鲁棒性更好。在随后的研究中我们希望将方法应用于流形优化，并且尽可能多地进行不同种类的数值实验来扩大我们算法的应用范围。

文章引用

张原浩. MCP正则优化问题的高效二阶算法
Efficient Second-Order Algorithm for MCP Regular Optimization Problems[J]. 应用数学进展, 2022, 11(10): 7173-7184. https://doi.org/10.12677/AAM.2022.1110761

参考文献