基于稀疏逻辑回归的加速分裂增广拉格朗日法

Accelerated Splitting Augmented Lagrange Method Based on Sparse Logistic Regression

摘要: 近几年来，稀疏优化在最小二乘、压缩感知尤其在统计学的逻辑回归模型中均得到了极大的研究及关注。逻辑回归是经典的分类方法，广泛应用于数据挖掘、机器学习、生物信息学。现研究带有l0范数约束的稀疏逻辑回归问题，为了求解此类问题，提出了加速分裂增广拉格朗日法(ASALM)。首先我们通过变量分裂，将稀疏逻辑回归问题转化为带有一个等式约束和一个基数约束的优化问题。接着在分裂增广拉格朗日法(SALM)中运用Nesterov加速技术，来提高算法的全局收敛性。然后在算法迭代中交替地求解一个无约束凸优化问题和一个带基数约束的二次优化问题，同时证明了ASALM的收敛性。最后，通过数值实验说明了ASALM的有效性和可行性，并且证明了ASALM的收敛速度比SALM更快。

Abstract: In recent years, sparse optimization has received great research and attention in the least squares, compressed sensing, especially in statistical logistic regression models. Logistic regression is a classic classification method, which is widely used in data mining, machine learning, and bioinformatics. Now we study the sparse logistic regression problem with l0 norm constraints. In order to solve such problems, the accelerated split augmentation Lagrangian method (ASALM) is proposed. First, we use variable splitting to transform the sparse logistic regression problem into an optimization problem with an equality constraint and a cardinality constraint. Next, Nesterov’s acceleration technology is used in the split augmented Lagrangian method (SALM) to improve the global convergence of the algorithm. Then, an unconstrained convex optimization problem and a quadratic optimization problem with cardinality constraints are alternately solved in the algorithm iteration; at the same time, the convergence of ASALM is proved. Finally, numerical experiments demonstrate the effectiveness and feasibility of ASALM, and prove that the convergence speed of ASALM is faster than SALM.

1. 引言

在如今的大数据时代下，稀疏优化的应用十分广泛，机器学习、压缩感知、分子生物学等领域都有很大的研究与发展。近几年来，稀疏优化在统计学的逻辑回归模型中得到了极大的研究及关注。逻辑回归分析是一种有效的数据处理分类方法，被广泛应用于文本分类 [1] [2]、机器学习、生物信息学 [3] [4] [5] [6] 等领域。经典的逻辑回归是一种线性的分类方法，具有训练速度快、模型的可解释性好等特点。本文考虑的是二项逻辑回归模型，通过给定的训练数据集，运用极大似然估计法估计模型参数得到逻辑回归模型。

在实际应用中，需要从大量特征中选择出对模型起关键作用的特征，排除与模型不相关或对模型作用很小的特征。同时也需要控制逻辑回归模型中参数的稀疏性，起到稀疏优化的效果。为了不造成特征泛滥，即需要控制模型中参数w非零元素的个数 ${\Vert w\Vert }_0$。而带有l₀范数的稀疏优化模型是NP-难问题，并且模型是非凸的，这给我们的计算带来了挑战，为了克服这个困难，一般的方法是加入正则项来进行松弛，将l₀范数松弛为l₁范数是常用的方法。带有l₁范数的正则化模型有着很好的凸性，能够防止过拟合，同时能够满足得到模型稀疏解，而被广泛应用在各个领域。例如带有l₁范数的线性回归模型通常又记为Lasso模型 [7]，在疾病分类中有着广泛的应用。添加l₁范数后的稀疏逻辑回归模型是一个凸优化问题，可以运用经典的凸优化方法去求解，例如一阶算法 [8]、迭代算法 [9]、内点算法 [10]。

在文献 [11] 中，通过对Nesterov加速 [12] 的研究，提出了加速增广拉格朗日法(AALM)，建立了该方法的迭代复杂度为 $O\left(1/k^2\right)$，并且又推导出了增广拉格朗日法(ALM)的迭代复杂度为 $O\left(1/k\right)$，证明了AALM比ALM具有更好的加速效果。本文将Nesterov加速方法应用到SALM框架中，提高了算法的全局收敛速度。针对ASALM框架中无约束凸优化子问题，运用拟牛顿法(BFGS算法)求解。针对ASALM框架中带基数约束的子问题，通过等价变换得到带l₀范数约束二次优化问题的显式解。

本文首先介绍了逻辑回归的基础知识，建立了逻辑回归模型和稀疏逻辑回归模型，并给出了最优性条件；然后针对稀疏逻辑回归模型，提出了加速分裂增广拉格朗日算法(ASALM)；又给出了该算法的一阶收敛性证明；然后，通过数值实验说明了ASALM的有效性和可行性。最后，对全文进行了总结。

2. 稀疏逻辑回归

稀疏逻辑回归基于逻辑回归模型，并在逻辑回归模型基础上添加l₀范数约束来控制模型中非零参数的个数。在实际应用中，需要从大量特征中选择出对模型起关键作用的主要特征。例如：在文字识别中，为了选择出对文字有鉴别意义的特征，需要在逻辑回归模型中添加l₀范数约束来控制模型中非零参数的个数。下面将对稀疏逻辑回归模型给出简单的回顾。

2.1. 二项逻辑回归模型

给定 $x={\left(x_1,\cdots ,x_n\right)}^{\text{T}}$ 表示输入变量，其中 $x_i$ 表示x的第i个属性(特征)， $y\in \left\{-1,+1\right\}$ 表示对应输出变量， $y_i=-1$ 表示第i个样本 $x^i$ 输出为负类， $y_i=+1$ 表示第i个样本 $x^i$ 输出为正类。对于给定的输入样本 $x\in {ℝ}^n$，输出标签为 $y\in \left\{-1,+1\right\}$ 的条件概率为：

$P\left(y=-1|x\right)=\frac{1}{1+{\text{e}}^{w^{\text{T}}x+b}},$

$P\left(y=\text{+}1|x\right)=\frac{{\text{e}}^{w^{\text{T}}x+b}}{1+{\text{e}}^{w^{\text{T}}x+b}},$

其中 $w\in {ℝ}^n,b\in ℝ$ 分别为权值向量与偏差，也为模型的待估计参数。给定训练数据集 $D=\left\{\left(x^1,y_1\right),\cdots ,\left(x^N,y_N\right)\right\}\subset {ℝ}^n\times \left\{-1,+1\right\}$，其中 $y_i=-1$ 表示第i个样本 $x^i$ 输出为负类， $y_i=+1$ 表示第i个样本 $x^i$ 输出为正类。采用极大似然估计法估计参数 $w,b$，等价变形后，通过求解极小化平均逻辑似然函数估计模型参数，得到如下逻辑回归模型。

P(1) $\min l_{avg}\left(w,b\right)=\min \frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\log }\left(1+\exp \left(-y_i\left(w^{\text{T}}{x}^i+b\right)\right)\right),$

其中 $w\in {ℝ}^n$， $b\in ℝ$ 为待估计参数。

2.2. 稀疏逻辑回归模型

稀疏逻辑回归是基于经典逻辑回归模型构建，并在经典逻辑回归模型基础上添加l₀范数约束来控制模型中非零参数的个数。因此稀疏逻辑回归模型构建为：

P(2) $\underset{w,b}{\min }{l}_{avg}\left(w,b\right),\text{s}\text{.t}\text{.}{\Vert w\Vert }_0\le K,$

其中 $\Vert \cdot \Vert$ 为欧式范数， $K>0$ 为控制w中非零元素个数的参数。带l₀范数约束的逻辑回归问题高度非线性，且是NP-难问题。通过添加l₁范数，稀疏逻辑回归模型转化为带l₁正则项的凸优化问题：

P(3) $\underset{w,b}{\min }{l}_{avg}\left(w,b\right)+\tau {\Vert w\Vert }_1,$

其中 $\tau >0$ 为正则化参数。凸优化问题可以用迭代算法 [9]、内点算法 [10]、一阶算法 [8] 等多种凸优化方法求解。

2.3. 最优性条件

下面给出稀疏逻辑回归的最优性条件 [13] [14]。令 $\alpha =\left(b,w\right)\in {ℝ}^{n+1}$，即 $l_{avg}\left(b,w\right)=l_{avg}\left(\alpha \right)$。若 ${\alpha }^{*}$ 是稀疏逻辑回模型的局部最优解，则存在 ${\delta }^{\text{*}}\in {ℝ}^{n+1}$ 与 ${\alpha }^{*}$ 满足如下的KKT条件：

$\begin{array}{l}\nabla l_{avg}\left({\alpha }^{*}\right)+{\delta }^{*}=0,\\ {\delta }^{*}=0,j\in J.\end{array}$ (1.1)

其中 $J^{*}=\left\{1\right\}\cup \left\{2\le j\le n+1:{\alpha }_j^{*}\ne 0\right\}$。

3. 加速分裂增广拉格朗日法(ASALM)

3.1. 变量分裂

稀疏逻辑回归模型是一个带有l₀范数约束的最优化问题，采用变量分裂的技术，将问题转化为带有一个等式约束和一个基数约束的最优化问题。由于w、b分别代表逻辑回归模型中权值向量和偏差。令 $\alpha =\left(b,w\right)$，则 $b={\alpha }_1,w={\alpha }_J$，其中 $J=\left\{2,\cdots ,n+1\right\}$。引入辅助变量 $\beta \in {ℝ}^n$，则稀疏逻辑回归模型等价转化为：

P(4) $\underset{\alpha ,\beta }{\min }{l}_{avg}\left({\alpha }_1,{\alpha }_J\right),\text{s}\text{.t}\text{.}{\alpha }_J=\beta ,{\Vert \beta \Vert }_0\le K.$

3.2. Nesterov加速

Nesterov提出一种用于梯度方法的加速策略，在Nesterov加速 [12] 算法中，引入一个或者多个辅助变量，并利用当前变量的梯度来更新参数。具体而言，在SALM框架中，拉格朗日乘子更新后得到 ${\hat{\gamma }}^k$，増加了对拉格朗日乘子预测并修正的步骤，如下：

$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2},$

${\gamma }^{k+1}={\hat{\gamma }}^k\text{+}\frac{t_k-1}{t_{k+1}}\left({\hat{\gamma }}^k-{\hat{\gamma }}^{k-1}\right)+\frac{t_k}{t_{k+1}}\left({\hat{\gamma }}^k-{\gamma }^k\right).$

3.3. 加速分裂增广拉格朗日法

关于等式约束 ${\alpha }_J=\beta$ 的增广拉格朗日函数形式为：

$L_{\rho }\left(\alpha ,\beta ,\gamma \right)=l_{avg}\left({\alpha }_1,{\alpha }_J\right)+{\gamma }^{\text{T}}\left(\beta -{\alpha }_J\right)+\frac{\rho }{2}{\Vert \beta -{\alpha }_J\Vert }^2,$

其中 $\gamma \in {ℝ}^n$ 为等式约束 ${\alpha }_J=\beta$ 对应的拉格朗日乘子， $\rho >0$ 为罚参数。亦即：

$L_{\rho }\left(\alpha ,\beta ,\gamma \right)=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\log \left(1+\exp \left(-B_i\alpha \right)\right)}+\frac{\rho }{2}{\Vert {\alpha }_J-\beta -\frac{1}{\rho }\cdot \gamma \Vert }^2-\frac{1}{2\rho }{\Vert \gamma \Vert }^2,$

其中 $B_i$ 为矩阵B的第i行，

$B=\left(\begin{array}{cc}{y}_1& y_1{\left(x^1\right)}^{\text{T}}\\ y_2& y_2{\left(x^2\right)}^{\text{T}}\\ ⋮& ⋮\\ y_N& y_N{\left(x^N\right)}^{\text{T}}\end{array}\right)\in {ℝ}^{N\times \left(n+1\right)}.$ (2.1)

下面给出加速分裂增广拉格朗日法(ASALM)算法步骤：

初始步：给定允许误差 $\epsilon \ge 0$，罚参数 $\rho >0$，初始化 $\hat{\gamma }\in {ℝ}^n$，初始化 $\beta$，令 ${\gamma }^1={\hat{\gamma }}^0$， $t_1=1$， $k=1$。

迭代步：判断精度，若 $\Vert {\alpha }_J^k-{\beta }^k\Vert \le \epsilon$ 输出 ${\alpha }^k$，停止计算；

令 $t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}$，计算：

${\alpha }^{k\text{+}1}=\arg \underset{\alpha \in {ℝ}^{n+1}}{\min }{L}_{\rho }\left(\alpha ,{\beta }^k,{\gamma }^k\right),$ (2.2)

${\beta }^{k\text{+}1}=\arg \underset{{\Vert \beta \Vert }_0\le K}{\min }{L}_{\rho }\left({\alpha }^{k\text{+}1},\beta ,{\gamma }^k\right),$ (2.3)

${\hat{\gamma }}^k={\gamma }^k+\rho \left({\beta }^{k\text{+}1}-{\alpha }_J^{k+1}\right),$ (2.4)

${\gamma }^{k+1}={\hat{\gamma }}^k\text{+}\frac{t_k-1}{t_{k+1}}\left({\hat{\gamma }}^k-{\hat{\gamma }}^{k-1}\right)+\frac{t_k}{t_{k+1}}\left({\hat{\gamma }}^k-{\gamma }^k\right),k:=k+1.$ (2.5)

由ASALM算法迭代步可知，子问题(2.2)是一个无约束光滑凸优化问题，可以用拟牛顿法求解，其核心思想是用不含二阶导数的矩阵 $U_k$ 替代牛顿法中的 $H_k^{-1}$ (Hesse矩阵的逆矩阵)，然后沿搜索方向 $-U_k\cdot g_k$ 做一维搜索，根据不同的 $U_k$ 构造方法有不同的拟牛顿法。

令 $f\left(\alpha \right)=L_{\rho }\left(\alpha ,{\beta }^k,{\gamma }^k\right)$，具体拟牛顿法(BFGS法)步骤如下：

初始步：给定初始点 ${\alpha }^{\left(1\right)}$，允许误差 $\epsilon$，设置 $B_1^{-1}$ (通常取单位阵)， $k=1$。

迭代步：1) 记 $g_k=\nabla f\left({\alpha }^{\left(k\right)}\right)$，计算搜索方向 $d^{\left(k\right)}=-B_k^{-1}\cdot g_k$，

2) 从点 ${\alpha }^{\left(k\right)}$ 出发，沿着 $d^{\left(k\right)}$ 做一维搜索，获得最优步长并更新参数：

$\begin{array}{l}{\lambda }_k=\arg \underset{\lambda }{\min }f\left({\alpha }^{\left(k\right)}+\lambda \cdot d^{\left(k\right)}\right),\\ {\alpha }^{\left(k+1\right)}={\alpha }^{\left(k\right)}+{\lambda }_k\cdot d^{\left(k\right)},\end{array}$

3) 判断精度，若 $\Vert g_{k+1}\Vert <\epsilon$，则停止迭代，否则转(4)；

4) 计算 $\Delta g=g_{k+1}-g_k,\Delta \alpha ={\alpha }^{\left(k+1\right)}-{\alpha }^{\left(k\right)}$。更新 $B^{-1}$ ：

$B_{k+1}^{-1}=\left(I_{n+1}-\frac{\Delta \alpha \Delta g^{\text{T}}}{\Delta {\alpha }^{\text{T}}\Delta g}\right)B_k^{-1}\left(I_{n+1}-\frac{\Delta g\Delta {\alpha }^{\text{T}}}{\Delta {\alpha }^{\text{T}}\Delta g}\right)+\frac{\Delta \alpha \Delta {\alpha }^{\text{T}}}{\Delta {\alpha }^{\text{T}}\Delta g},$

5) $k=k+1$，转1)。

子问题(2.3)是一个带有基数约束的最优化问题，通过等价变换可以转化为带有基数约束的二次优化问题如下：

$\underset{\beta }{\min }{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left({\beta }_i-{\xi }_i\right)}^2},$

$\text{s}\text{.t}\text{.}{\Vert \beta \Vert }_0\le K.$

其中 $\xi ={\alpha }_J-\gamma /\rho$，且问题的最优解 [14] 为如下形式： $y_i^{\ast }={\xi }_i$， $i\in {\rm M}$。 $y_i^{\ast }=0$， $i\notin {\rm M}$。其中 ${\rm M}\subset \left\{1,2,3,\cdots ,n\right\}$ 表示向量 ${\xi }^2$ 的前K个最大元素所对应的指标集。

3.4. ASALM算法的收敛性

ASALM算法求解稀疏逻辑回归问题的一阶收敛性 [13] 如下：

令 $\left(\hat{\alpha },\hat{\beta },\hat{\gamma }\right)$ 为由ASALM算法产生点列 ${\left\{\left({\alpha }^k,{\beta }^k,{\gamma }^k\right)\right\}}_{k=1}^{\infty }$ 的聚点。假设 ${\left\{{\gamma }^k\right\}}_{k=1}^{\infty }$ 有界，且

$\underset{k\to \infty }{\lim }\left[\left({\alpha }^{k+1},{\beta }^{k+1},{\gamma }^{k+1}\right)-\left({\alpha }^k,{\beta }^k,{\gamma }^k\right)\right]=0.$

则 ${\hat{\alpha }}^k$ 满足问题P(4)的一阶最优性条件(1.1)。

4. 数值实验

本章节主要通过数值实验来验证加速分裂增广拉格朗日算法(ASALM)的可行性和有效性。我们将所要研究的稀疏逻辑回归问题应用于随机合成的模拟数据，通过对实验结果的数据分析和绘图分析客观地描述了加速分裂增广拉格朗日算法的性能与效率。本次实验的代码由MATLAB R2019a编写，并且在Intel(R) Core(TM) i5-9300H 2.4GHz CPU, 8GB RAM的个人笔记本电脑上进行了数值实验。

我们考虑下列稀疏逻辑回归问题：

$\underset{\alpha ,\beta }{\min }\left\{\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\log \left(1+\exp \left(-B_i\alpha \right)\right)}\right\},\text{s}\text{.t}\text{.}{\alpha }_J=\beta ,{\Vert \beta \Vert }_0\le K.$ (3.1)

其中 $\alpha =\left({\alpha }_1,{\alpha }_J\right)$，B是一个 $N\times \left(n+1\right)$ 矩阵如上(2.1)所示。设置罚参数 $\rho =0.5$， $k=1$， $t=1$，初始化 ${\beta }^1,{\hat{\gamma }}^0$ 为零向量，即 ${\gamma }^1$ 亦为零向量。重新定义矩阵B的维度尺寸为 $n\times m$，定义 $r^{k+1}={\alpha }_J^{\text{k}+1}-{\beta }^{k+1}$ 和 $s^{k+1}=-\rho \left({\beta }^{k+1}-{\beta }^k\right)$ 分别为原始残差和对偶残差，设置绝对误差 ${\epsilon }^{abs}={10}^{-6}$ 和相对误差 ${\epsilon }^{rel}={10}^{-3}$，实验中的迭代停止准则为 $\Vert r^k\Vert \le {\epsilon }^{pri}$ 和 $\Vert s^k\Vert \le {\epsilon }^{dual}$，其中 ${\epsilon }^{pri}=\sqrt{n}{\epsilon }^{abs}+{\epsilon }^{rel}\max \left\{\Vert {\alpha }^k\Vert ,\Vert {\beta }^k\Vert \right\}$， ${\epsilon }^{dul}=\sqrt{n}{\epsilon }^{abs}+{\epsilon }^{rel}\Vert {\gamma }^k\Vert$。

为了直观地描述ASALM的性能与效率，我们将ASALM与分裂增广拉格朗日算法(SALM)进行了实验对比，测试了四种不同维度的应用实例，即 $\left(n,m\right)=\left(500,500\right),\left(500,1000\right),\left(1000,1000\right),\left(1000,1500\right)$，数值实验结果见表1。在四种不同的维度条件下，表1中分别比较了两种算法的迭代次数、CPU运算时间和目标函数值。从表1中可以分析得出ASALM比SALM有着更少的迭代次数和CPU运算时间。因此，针对问题(3.1) ASALM是可行的且比SALM更有效。

为了进一步观察ASALM和SALM的收敛性，我们可视化了两种算法收敛性的变化，见如下图1和图2。我们观察了 $\left(n,m\right)=\left(500,500\right),\left(500,1000\right)$ 两种条件下原始残差和目标函数值随着迭代次数(k)的变化情况，图1和图2都展示了ASALM的原始残差比SALM下降的更快，且ASALM比SALM更快达到最优目标函数值。因此，ASALM的收敛速度比SALM更快，且性能与效率更好。

Table 1. Comparison of two algorithms based on simulated data

表1. 基于模拟数据的两种算法比较

Figure 1. n = 500, m = 500

图1. n = 500, m = 500

Figure 2. n = 500, m = 1000

图2. n = 500, m = 1000

5. 结论与展望

本文针对带有l₀范数约束的逻辑回归问题，提出了加速分裂增广拉格朗日法求解此类问题。该算法在迭代中交替求解两个子问题，借助BFGS算法求解其中一个子问题，并通过等价变形求解出另一个子问题的显式解。在SALM框架中，运用Nesterov加速技术，提高算法的全局收敛性，给出了ASALM算法的一阶收敛性证明。最后，通过数值实验说明了ASALM的收敛速度比SALM更快，且性能与效率更好。此类方法亦可以应用到带有范数约束的多项逻辑回归问题上，在提高全局收敛性的性能上优于其它凸优化方法。