Stiefel流形上非光滑优化的一种带外推的可变度量邻近梯度算法

A Variable Metric Proximal Gradient Method with Extrapolation for Nonsmooth Optimization over the Stiefel Manifold

摘要: 本文针对Stiefel流形上一类目标函数为光滑损失函数与非光滑函数之和的非凸优化问题，提出了一种基于收缩的可变度量惯性邻近梯度算法。所提出的算法在已有的加速黎曼邻近梯度算法基础上，引入了对角Barzilai-Borwein类步长策略，该策略可以更好的捕获问题的局部几何信息，进一步加速算法的收敛。理论上，证明了算法全局收敛到稳定点。最后，本文给出了稀疏主成分分析问题的数值结果，验证了该方法的有效性。

Abstract: In this paper, we propose a retraction based variable metric inertial proximal gradient method for solving a class nonconvex optimization problem over the Stiefel manifold whose objective function is the summation of a smooth cost function and a nonsmooth function. Based on existing inertial Riemannian proximal gradient method, the proposed method introduces a metric changing called diagonal Barzilai-Borwein step-size strategy at each iteration, which can better capture the local geometric of this class problem and accelerate the convergence of the algorithm. Theoretically, we show that the proposed method globally converges to a stationary point. Numerical results on solving sparse PCA problem is reported to demonstrate the efficiency of the proposed method.

1. 引言

本文主要考虑以下具有Stiefel流形约束的一类非凸非光滑优化问题：

$\begin{array}{l}\min F\left(x\right):=f\left(x\right)+g\left(x\right)\\ \text{s}\text{.t}\text{.}x\in \mathcal{M}\end{array}$ (1)

其中函数 $f:\mathcal{M}\to ℝ$ 是可微的但可能是非凸的，函数 $g:\mathcal{M}\to ℝ$ 是连续的但也可能是非光滑的。可行集 $\mathcal{M}:={\mathcal{S}}_{n,p}=\left\{X:X\in {ℝ}^{n\times p},X^{\text{T}}X=I_p\right\}$ 称为Stiefel流形，也称为正交约束，其中 $I_p$ 记为p阶单位阵( $p\le n$)。

Stiefel流形上的非光滑优化问题由于其在各个领域的广泛应用引起了许多研究者的注意。例如，在统计与数据科学中的稀疏主成分分析问题 [1]；物理中的压缩模问题 [2]；图像领域中的盲去卷积问题 [3] 以及机器学习中的无监督特征选取 [4] 等问题，其本质都是Stiefel流形上的非光滑优化问题。更多的相关应用，我们推荐读者参考文献 [5]。

通常情况下，由于其目标函数的不可微性和流形约束的非凸性，问题(1)是较难求解的。对于部分问题已经证明是NP-难的 [6]。求解问题(1)的数值算法相较于目标函数是光滑形式的算法是有限的，目前已有的算法包括次梯度类算法 [7] [8]、算子分裂算法 [9] [10] 以及邻近梯度算法 [11] [12]。然而，上述所提到的算法存在着寻找下降方向成本高、缺乏收敛性分析、对参数敏感等问题。众所周知，对于目标函数为光滑损失函数与非光滑连续函数之和的复合结构，当其对应的邻近算子容易求解时，邻近梯度算法是求解该类问题的最有效方法之一。然而，由于问题(1)中Stiefel流形的约束，此时邻近算子是没有显式解的。针对这一问题，最近Chen等 [13] 提出了一种基于收缩的邻近梯度算法，其中下降方向是由限制在Stiefel流形切空间中的邻近子问题决定。尽管该子问题没有显示解，作者通过将其转换为一个非线性方程系统，进而采用正则的半光滑牛顿法 [14] 高效求解。随后，Huang等 [15] 将快速迭代收缩阈值算法(FISTA)从欧氏空间推广到黎曼空间中。进一步，Huang等 [16] 在2021年针对一般的黎曼流形研究了邻近梯度算法及其Nesterov加速版本，作者进一步借助Kurdyka-Łojasiewicz不等式性质分析了所提出算法的局部收敛率。虽然上述文献中均假设目标函数中光滑部分是梯度Lipschitz连续的，然而在实际问题中，其Lipschitz常数的精确计算并不是容易的。

受文献 [17] [18] 的启发，本文基于上述研究成果，提出了一种可变度量的惯性邻近梯度算法求解问题(1)，所提出的算法采用对角化的Barzilai-Borwein (BB)步长估计Lipschitz常数，同时结合Nseterov动量项进一步加速算法收敛，并给出收敛性分析。最后，在稀疏主成分分析问题中的数值结果验证了我们的方法是有效的。

2. 预备知识

本节引入一些流形优化的基本定义和概念。

定义 1 (收缩映射) 一个光滑映射 ${\mathrm{Retr}}_X:T_X\mathcal{M}\to \mathcal{M}$ 满足

1) ${\mathrm{Retr}}_X\left(0\right)=X,\forall X\in \mathcal{M}$，其中0定义为 $T_X\mathcal{M}$ 的零元素；

2) 对于任意的 $X\in \mathcal{M}$，有如下极限成立：

$\underset{0\leftarrow \xi \in T_X\mathcal{M}}{\lim }\frac{{\Vert {\mathrm{Retr}}_X\left(\xi \right)-\left(X+\xi \right)\Vert }_F}{{\Vert \xi \Vert }_F}=0$,

则 ${\mathrm{Retr}}_X$ 被称为收缩映射。

定义2 (广义Clarke次微分)对于一个定义在 $\mathcal{M}$ 上的局部Lipschitz函数F，则F在 $X\in \mathcal{M}$ 处方向V上的黎曼方向导数定义为：

$F^{\circ }\left(X,V\right)=\underset{Y\to X,t\downarrow 0}{\lim }\sup \frac{F\circ {\varphi }^{-1}\left(\varphi \left(Y\right)+tD\varphi \left(X\right)\left[V\right]\right)-F{\varphi }^{-1}\left(\varphi \left(Y\right)\right)}{t}$,

其中 $\left(\varphi ,U\right)$ 是在X处的坐标卡。则函数F在 $X\in \mathcal{M}$ 处的广义梯度或者Clarke次微分表示为 $\hat{\partial }F\left(X\right)=\left\{\xi \in T_X\mathcal{M}:\langle \xi ,V\rangle \le F^{\circ }\left(X,V\right)\forall V\in T_X\mathcal{M}\right\}$

由于F是正则函数，本文有 $\hat{\partial }F\left(X\right)={\mathrm{Proj}}_{T_X\mathcal{M}}\left(\partial F\left(X\right)\right)$，其中 ${\mathrm{Proj}}_{T_X\mathcal{M}}$ 是正交投影。进一步我们可以得出问题(1)的一阶最优性条件，即 $0\in gradf\left(X\right)+{\mathrm{Proj}}_{T_X\mathcal{M}}\left(\partial g\left(X\right)\right)$，其中 $\partial (·)$ 表示为欧式空间中的次微分。

3. 可变度量的惯性邻近梯度算法

本节中，我们提出Stiefel流形约束优化问题可变度量的惯性邻近梯度算法。

为了满足约束，本文的算法通过更新迭代

$X_{k+1}={\mathrm{Retr}}_{X_k}\left({\alpha }_k{d}_k\right)$ (2)

来保持可行性，其中 ${\alpha }_k$ 是由Armijo条件决定的， $d_k$ 是搜索方向。如同文献 [13] 中所示，搜索方向 $d_k$ 是由如下限制在切空间中的邻近子问题所决定

$d_k=\underset{d\in T_{X_k}\mathcal{M}}{\arg \min }\langle gradf\left(X_k\right),d\rangle +\frac{L}{2}{\Vert d\Vert }_F^2+g\left(X_k+d\right)$, (3)

其中L为可微函数f的梯度Lipschitz常数。

具体地，本文类似于无约束优化中的情况，首先给出如下的长和短BB步长：

${\alpha }_{k+1}^{BB1}=\underset{\alpha \in ℝ}{\arg \min }{\Vert {\alpha }^{-1}{S}_k-Y_k\Vert }_F=\frac{tr\left(S_k^{\text{T}}{S}_k\right)}{|tr\left(S_k^{\text{T}}{Y}_k\right)|}$,

和

${\alpha }_{k+1}^{BB2}=\underset{\alpha \in ℝ}{\arg \min }{\Vert S_k-\alpha Y_k\Vert }_F=\frac{tr\left(S_k^{\text{T}}{Y}_k\right)}{|tr\left(Y_k^{\text{T}}{Y}_k\right)|}$

其中 $S_k=X_{k+1}-X_k$， $Y_k=d_{k+1}-d_k$。这里下降方向 $d_{k+1}$ 和 $d_k$ 分别是限制在不同切空间中邻近子问题的解。为了更好的获取函数f的Hessian信息，采用如下在第k次迭代中计算得出的度量 $U^k=Diag\left(u^k\right)$， $u^k=\left[u_1^k,\cdots ,u_n^k\right]\in {ℝ}^n$，其中

$u_i^k=\{\begin{array}{l}\frac{1}{{\alpha }_k^{BB1}}\frac{S_k^i{Y}_k^i+\upsilon u_i^{k-1}}{{\left(S_k^i\right)}^2+\upsilon }<\frac{1}{{\alpha }_k^{BB1}}\\ \frac{1}{{\alpha }_k^{BB2}}\frac{S_k^i{Y}_k^i+\upsilon u_i^{k-1}}{{\left(S_k^i\right)}^2+\upsilon }>\frac{1}{{\alpha }_k^{BB2}}\\ \frac{S_k^i{Y}_k^i+\upsilon u_i^{k-1}}{{\left(S_k^i\right)}^2+\upsilon }否则\end{array}$ (4)

此时，邻近子问题(3)就变为

$d_k=\underset{d\in T_{X_k}\mathcal{M}}{\arg \min }\langle gradf\left(X_k\right),d\rangle +\frac{L}{2}\langle d,U^{k}d\rangle +g\left(X_k+d\right)$ (5)

该子问题仍采用半光滑牛顿法求解。接下来，我们给出求解问题(1)的可变度量惯性邻近梯度算法。

算法1. 可变度量惯性邻近梯度法

VM-AManPG算法：

步0：初始化 $X_0\in \mathcal{M}$，线搜索参数 $\delta$ 、 $\sigma \in \left(0,1\right)$，超参数 $\upsilon >0$，保护步中正整数M， $t_0=1$。

步1： $Y_0=X_0$， $Z_0=X_0$ ；

步2：for $k=0,\cdots$ do

步3： ifmod(k, M)=0 then

步4： 调用算法2： $\left[Z_{k+M},X_k,Y_k,t_k\right]=\mathrm{Alg}2\left(Z_{k+M},X_k,Y_k,t_k,F\left(X_k\right)\right)$ ；

步5： end if

步6：通过半光滑牛顿法求解子问题(5)得出下降方向 $d_k$ ；

步7：步通过式(2)更新 $X_{k+1}$ ；

步8： $t_{k+1}=\frac{\sqrt{4t_k^2+1}+1}{2}$ ；

步9：计算 $Y_{k+1}={\mathrm{Retr}}_{X_{k+1}}\left(\frac{1-t_k}{t_{k+1}}{\mathrm{Retr}}_{X_{k+1}}^{-1}\left(X_k\right)\right)$ ；

步10：end for

注1. 其中VM-AManPG算法步骤9中的收缩映射 ${\mathrm{Retr}}_{(·)}$ 本文考虑使用极分解，

${\mathrm{Retr}}_X^{polar}\left(\xi \right)=\left(X+\xi \right){\left(I_p+{\xi }^{\text{T}}\xi \right)}^{-1/2}$，且该收缩的逆是存在的。众所周知，惯性邻近梯度法是非单调的，

基于文献 [15]，本文同样采用重启策略使得算法VM-AManPG是满足下降性的。

算法2. 算法1.的保护步

VM-AManPG算法：

步0：输入 $\left(Z_k,X_k,Y_k,t_k,F\left(X_k\right)\right)$ ；

步1：计算子问题(5)得到下降方向 $d_{Z_k}$ ；

步2：令 $\alpha =1$ ；

步3： while $F\left({\mathrm{Retr}}_{Z_k}\left(\alpha d_{Z_k}\right)\right)>F\left(Z_k\right)-\sigma \alpha {\Vert d_{Z_k}\Vert }_F^2$ do

步4： $\alpha =\delta \alpha$ ；

步5： end while

步6：若 $F\left({\mathrm{Retr}}_{Z_k}\left(\alpha d_{Z_k}\right)\right)<F\left(X_k\right)$ 则

步7： $X_k={\mathrm{Retr}}_{Z_k}\left(\alpha d_{Z_k}\right)$ ； $Y_k={\mathrm{Retr}}_{Z_k}\left(\alpha d_{Z_k}\right)$， $t_k=1$ ；

步8：否则 $X_k$， $Y_k$， $t_k$ 若保持不变；

步9：输出 $Z_{k+M}=X_k$。

4. 收敛性分析

本节中分析了VM-AManPG算法的收敛性质。在正式给出结论之前，本文需要做出如下的假设。

假设1. 函数f是可微的，且其梯度 $\nabla f$ 是Lipschitz连续的有常数L。函数g是凸非光滑的，并且是Lipschitz连续的。

假设2. 函数F是强制的，即当 ${\Vert X\Vert }_F\to \infty$ 时 $F\left(X\right)\to +\infty$。

假设3. 存在两个正常数 $0<\gamma <\bar{\gamma }$，使得对角矩阵U在任意点 $X_k$ 处的特征值满足在 $\gamma ,\bar{\gamma }$ 范围内。

引理1：序列 $\{Z_k\},\{d_k\}$ 分别是由算法VM-AManPG生成的迭代序列，当假设1成立时，存在常数 $\bar{\alpha }$ 使得对任意的 $0<\alpha \le \min \left(1,\bar{\alpha }\right)$，满足

$F\left({\mathrm{Retr}}_{Z_k}\left(\alpha d_k\right)\right)-F\left(Z_k\right)\le -\frac{\alpha L}{2}{\Vert d_k\Vert }_F^2$

上述引理说明算法1是有定义的。事实上，算法1中序列 $\{Z_k\}$ 的下标为 $k+M$。当本文中所采用的对角BB步长矩阵 $U^k$ 为单位阵时，子问题(5)退化为文献 [13] 中所求解的子问题，此时，文献 [13] 中的结论在本文中都成立。进一步假设3成立的条件下，无论 $U^k$ 是否为单位阵，上述引理都是成立的。

定理1：在假设1、2和假设3成立的条件下，序列 $\{Z_k\}$ 是由算法1产生的迭代序列。令 $Z_{*}$ 是序列 $\{Z_k\}$ 的任意聚点，则有

$0\in {\mathrm{Proj}}_{T_Z\mathcal{M}}\partial F(\; Z\; *\; )$

即序列 $\{Z_k\}$ 的任意聚点是问题(1)的一个稳定点。

证明：在假设1和假设2条件下，由于Stiefel流形是紧集，则次水平集 ${\Omega }_{X_0}=\{X\in \mathcal{M}|F\left(X\right)\le F\left(X_0\right)\}$

是紧的，由于假设函数F是强制的，所以次水平集 ${\Omega }_{X_0}$ 是有界的。进一步由于F是连续函数，则F在次水平集 ${\Omega }_{X_0}$ 是有界的。

当下降方向 $d_k=0$ 时，根据子问题(5)的最优性条件有

$0\in gradf\left(Z_k\right)+\frac{1}{t}{U}^k{d}_k+{\mathrm{Proj}}_{T_Z\mathcal{M}}\partial g\left(Z_k+d_k\right)$ (6)

成立，这恰好是原问题(1)的一阶必要性条件。接下来我们就说明序列 $\{d_k\}$ 是满足 $\underset{k\to \infty }{\lim }{\Vert d_k\Vert }_F=0$。

根据算法2中步骤3的Armijo条件，有 $F\left(Z_{k+1}\right)<F\left(Z_k\right)$，即函数序列 $\{F\left(Z_k\right)\}$ 是单调下降的，进一步，由于函数F是有下界的，则序列 $\{F\left(Z_k\right)\}$ 的极限存在，有

$\underset{k\to \infty }{\lim }F\left(Z_k\right)-F\left({\mathrm{Retr}}_{Z_k}\left(\alpha d_k\right)\right)=0$

结合引理1，则有 $\underset{k\to \infty }{\lim }{\Vert d_k\Vert }_F=0$ 成立。

由子问题(5)的最优性条件(6)，在式(6)两边分别加上 $gradf\left(Z_k+d_k\right)$，可得

$\begin{array}{l}gradf\left(Z_k+d_k\right)-gradf\left(Z_k\right)-\frac{1}{t}{U}^k{d}_k\\ \in {\mathrm{Proj}}_{T_Z\mathcal{M}}\partial g\left(Z_k+d_k\right)+gradf\left(Z_k+d_k\right)={\mathrm{Proj}}_{T_Z\mathcal{M}}\partial F\left(Z_k+d_k\right)\end{array}$ (7)

且有

$\begin{array}{l}{\Vert gradf\left(Z_k\right)-gradf\left(Z_k+d_k\right)+\frac{1}{t}{U}^k{d}_k\Vert }_F\\ ={\Vert {\mathrm{Proj}}_{T_Z\mathcal{M}}\nabla f\left(Z_k\right)-{\mathrm{Proj}}_{T_Z\mathcal{M}}\nabla f\left(Z_k+d_k\right)+\frac{1}{t}{U}^k{d}_k\Vert }_F\\ \le \left(L+\frac{1}{t}\right){\Vert d_k\Vert }_{U^k}\to 0\end{array}$ (8)

当 $k\to \infty$ 时，其中 ${\mathrm{Proj}}_{T_Z\mathcal{M}}(·)$ 为光滑映射，满足范数不等式。

令 $\{Z_{k_i}\}$ 为 $\{Z_k\}$ 的子序列且收敛到 $Z_{*}$。接下来我们证明 $0\in {\mathrm{Proj}}_{T_Z\mathcal{M}}\partial F\left(Z_{*}\right)$。根据式(7)，存在序列 $\{{\tau }_k\}\in N_{Z_k}\mathcal{M}$，其中 $N_{Z_k}\mathcal{M}$ 为流形 $\mathcal{M}$ 在点 $Z_k$ 处的法空间，于是有

$gradf\left(Z_k+d_k\right)-gradf\left(Z_k\right)-\frac{1}{t}{U}^k{d}_k+{\tau }_k\in \partial F\left(Z_k+d_k\right)$

进一步结合当 $i\to \infty$ 时， $\{Z_{k_i}\}\to Z_{*}$，则有

$gradf\left(Z_{k_i}+d_{k_i}\right)-gradf\left(Z_{k_i}\right)-\frac{1}{t}{U}^{k_i}{d}_{k_i}+{\tau }_{k_i}\in \partial F\left(Z_{k_i}+d_{k_i}\right)$

由于目标函数F在紧集 ${\Omega }_{X_0}$ 上是连续的，则存在常数 $G>0$ 使得 ${\max }_{Z\in {\Omega }_{X_0}}{\max }_{\omega \in \partial F\left(Z\right)}{\Vert \omega \Vert }_F<G$，于是有 ${\Vert {\tau }_{k_i}\Vert }_F<G$，存在收敛子列 $\{{\tau }_{k_{i_j}}\}$ 且极限点为 ${\tau }_{*}$。则当 $j\to \infty$ 时，有

$gradf\left(Z_{k_{i_j}}+d_{k_{i_j}}\right)-gradf\left(Z_{k_{i_j}}\right)-\frac{1}{t}{U}^{k_{i_j}}{d}_{k_{i_j}}+{\tau }_{k_{i_j}}\to {\tau }_{*}$

以及 $Z_{k_{i_j}}+d_{k_{i_j}}\to Z_{*}$。因此有 $0\in \partial F\left(Z_{*}\right)$，进一步根据投影 ${\mathrm{Proj}}_{N_Z\mathcal{M}}$ 是光滑的，有

${\mathrm{Proj}}_{N_{Z_{k_{i_j}}}\mathcal{M}}{\tau }_{k_{i_j}}\to {\mathrm{Proj}}_{N_{Z_{{}_{*}}}\mathcal{M}}{\tau }_{*}$

所以，我们可以得出 $0\in {\mathrm{Proj}}_{T_{Z_{{}_{*}}}}\partial F\left(Z_{*}\right)$。

5. 数值实验

本节中，我们将在数据分析领域中的稀疏主成分分析问题上验证算法的有效性。所有的数值实验均是通过运行环境为64bit Ubuntu platform，CPU (Intel Core i5-5200U)为2.20GHz的Matlab R2019a实现。

稀疏主成分分析(SPCA)已经成为一种强大的数据分析技术，通过识别数据中的本地空间结构和消除不同时间尺度之间的歧义，提供对低秩结构的改进描述。该问题在文献 [1] 中经过离散化，可以转化为如下的模型：

$\begin{array}{l}\min -\mathrm{Tr}\left(X^{\text{T}}{A}^{\text{T}}AX\right)+\mu {\Vert X\Vert }_1\\ \text{s}\text{.t}\text{.}X\in S_{n,p},\end{array}$

其中矩阵 $A\in {ℝ}^{m\times n}$， ${\Vert X\Vert }_1={\displaystyle {\sum }_{ij}\left|X_{ij}\right|}$， $\mu >0$ 是正则化参数，可以控制上述模型的稀疏性。

本文通过与最新的算法ManPG-Ad [13]，AManPG [15]，ARPG [16] 进行比较来验证所提出算法的有效性。在数值实验中，算法VM-AManPG的参数设置如下：上述算法的终止条件均设置为

${\Vert d_k/t\Vert }_F^2\le \epsilon :={10}^{-8}np$。超参数 $\upsilon =2$，线搜索参数 $\delta =0.5$， $\sigma =0.001$，常数 $L=4{\Vert A\Vert }_F^2$，参数 $M=5$。另

外，算法ManPG-Ada，AManPG ，ARPG 的相关参数均保持为文献中的默认值。

在数值试验中，通过两类不同的数据测试算法。第一类是采取随机数据，即矩阵A和变量X都是随机生成的，对于数据矩阵A，其生成方式为：首先产生随机矩阵 $A=\mathrm{randn}\left(m,n\right)$，且使该矩阵的列均值为0，最后使其列的欧氏范数等于1。在所有的测试中均令 $m=50$。第二类是采用文献 [15] 中真实的DNA数据集，其变量维数 $n=24589$， $m=113$， $p=4$。

本文首先基于随机生成的数据，通过选取不同的变量维数n，p和正则化参数 $\mu$ 来测试算法VM-AManPG的性能。分别从表1、表2、表3中我们可以发现，VM-AManPG无论在CPU时间上还是在迭代次数上都明显的优于算法ManPG-Ada、AManPG和ARPG。数值结果表明，VM-AManPG具有明显的加速效果，并且说明了对角BB步长策略的有效性。

Table 1. Comparison on SPCA model with varies n with p = 5 and regular parameter μ = 0.5

表1. 在SPCA模型上不同变量维数n，其中维数 $p=5$，正则化参数 $\mu =0.5$ 地比较

Table 2. Comparison on SPCA model with varies p with n = 5000 and regular parameter μ = 0.5

表2. 在SPCA模型上不同变量维数p，其中维数 $n=5000$，正则化参数 $\mu =0.5$ 地比较

Table 3. Comparison on SPCA model with varies regular parameter μ with n = 3000 and p = 5

表3. 在SPCA模型上不同正则化参数 $\mu$，其中维数 $n=3000$， $p=5$ 地比较

接下来我们采用真实DNA数据集，通过选取不同的正则化参数 $\mu$ 评估算法性能。数值结果如表4所示，我们所提出的算法与已有的算法相比较，满足相同终止条件所需的CPU时间和迭代次数都具有明显的优势。

Table 4. Comparison on SPCA model with varies regular parameter μ , DNA data n = 24589 , m = 113 , p = 4

表4. 在SPCA模型上不同正则化参数 $\mu$，DNA 数据 $n=24589$， $m=113$， $p=4$ 地比较

6. 总结

本文针对一类在Stiefel流形上的非凸非光滑优化问题，提出了一种可变度量的惯性邻近梯度算法。通过将欧氏空间中的对角Barzilai-Borwein类步长策略推广到Stiefel流形上可以进一步加速算法收敛，在一定的假设条件下，可以证明算法全局收敛到原问题的一个稳定点。数值实验表明，无论是在随机数据还是在真实数据集中的数值结果，我们所提出的算法在求解效率上都具有一定的优势。