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

收稿日期：2020年9月1日；录用日期：2020年9月19日；发布日期：2020年9月27日

摘要

本文研究了分裂可行问题的条件梯度算法，该算法将求解迭代方向转化成求解一线性子问题，并以线搜索得到的步长作为凸因子，当前方向与上一步迭代点的凸组合作为新的迭代点。算法在迭代的更新步中不使用投影，并且得到的解有较好的稀疏性和低秩性。我们获得了算法的收敛性并给出数值实验对比分析了本文的算法与相关算法在同一算例下的表现情况，得到了良好的结果。

关键词

非线性分裂可行问题，条件梯度法，替代函数，稀疏约束集，IMRT问题

Solving Nonlinear Split Feasibility Problem via Conditional Gradient Method

Zhensheng Yu, Zilun Wang

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

Received: Sep. 1^st, 2020; accepted: Sep. 19^th, 2020; published: Sep. 27^th, 2020

ABSTRACT

In this paper, we consider a conditional gradient algorithm for the split feasibility problem. We get the iterative direction by solving a linear subproblem, use the step size obtained by the line search as the convex factor, and get the new iteration point by the convex combination of the current direction and the previous step iteration point. The algorithm does not use projection in the update step. And the obtained solution has good properties of sparsity and low rank. The numerical experiment compares the performance of the proposed algorithm with other algorithms under the same calculation example.

Keywords:Nonlinear Split Feasibility, Conditional Gradient Method, Surrogate Function, Sparse Constraint, IMRT

Copyright © 2020 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. 引言

本文主要研究欧几里得空间中的单集分裂可行问题(Nonlinear Split Feasibility Problem)，即找一个 $x^{*}\in R^n$ 使其满足

$x^{*}\in C$ 并且 $A{x}^{*}\in Q$ (1.1)

其中 $C\subseteq R_n$ 与 $Q\subseteq R_m$ 都为非空闭凸集，A代表矩阵或者线性变换。该问题由Censor等人 [1] 提出，它在放射治疗(IMRT)计划领域有着重要的作用 [2] [3]，在工程技术领域，例如信号和图像恢复，基因调控网络推论以及用于处理高维稀疏数据的LASSO和Dantzig Selector方法 [4] [5] [6] 中也有广泛应用。

循环投影法是最先被提出用于求解分裂可行问题的一种迭代算法，其迭代公式为：

$x^{k+1}=A^{-1}{\left(I+A{A}^{\text{T}}\right)}^{-1}\left(A{P}_C\left(x^k\right)+A{A}^{\text{T}}{P}_Q\left(A{x}^k\right)\right)$ (1.2)

这里， $P_C$ 和 $P_Q$ 分别表示向非空闭凸集C和Q上的正交投影。该方法显然不是很理想，因为在每次迭代中都需要计算矩阵的逆。当将该方法应用于大规模的实际问题时，求解逆矩阵非常困难，这也可能导致算法花费太长时间。因此，Byrne给出了著名的CQ算法 [7]。在之后的分裂可行问题算法的发展中，CQ算法得到更多的关注。假设C和Q分别是欧几里得或希尔伯特空间 $H_1$ 和 $H_2$ 上的闭合凸子集，并且A是有界线性算子。CQ算法的迭代格式如下

$x^{k+1}=P_C\left(x^k-\theta A^{*}\left(I-P_Q\right)A{x}^k\right)$ (1.3)

学者们还提出了一些其他的方法 [8] [9] [10] [11]，将分裂可行问题等价于不动点问题或变分不等式问题进行求解，因此应用于其等效问题的算法可以作适应性调整并应用于分裂可行问题，例如Krasnosel’skii-Mann (KM)算法 [12] [13] 在解决分裂可行问题方面取得了重大突破。Zhao将KM迭代形式应用于扰动CQ算法。由于KM迭代算法具有较弱的收敛性，因此作为一种特殊形式，CQ算法也具有较弱的收敛性。Xu [14] 改进了上述算法，并很好地解决了这一缺点，改进算法具有较强的收敛性。

考虑对分裂可行问题的一种情况：在欧式空间中存在集合 $C_i\subseteq R_n$，$Q_j\subseteq R_m$ 其中 $i=1,2,\cdots ,p$，$j=1,2,\cdots ,r$ 并满足(1.1)式条件，这个问题称为多集分裂可行性问题(MSSFP)，更进一步，我们设矢量值函数或非线性函数 $F:R_n\to R_m$ 代替A。本文中我们考虑双集合非线性分裂可行性问题，即

$x^{*}\in C\text{and}F\left(x^{*}\right)\in Q$ (1.4)

在2005年，Censor等人 [2] 提出了一种投影梯度算法，最小化一个使用最大最小法(Majorize-Minimize Method)构建的目标函数如下

$f\left(x\right)=\frac{1}{2}{\displaystyle \underset{i}{\sum }{v}_{i}dist{\left(x,C_i\right)}^2}+\frac{1}{2}{\displaystyle \underset{j}{\sum }{w}_{j}dist{\left[H\left(x\right),Q_j\right]}^2}$ (1.5)

其中 $dist\left(x,y\right)$ 表示向量x和y之间的欧氏距离。权重系数或者凸系数 $v_i$ 和 $w_j$ 为正数且满足 ${\displaystyle \sum v_i}+{\displaystyle \sum w_j}=1$。算法的迭代公式为：

$x^{k+1}=P_{\Omega }\left\{x^k+s\left({\displaystyle \underset{i=1}{\overset{t}{\sum }}{\alpha }_i\left(P_{C_i}\left(x^k\right)-x^k\right)}+{\displaystyle \underset{j=1}{\overset{r}{\sum }}{\beta }_j{A}^{\text{T}}\left(P_{Q_j}\left(A{x}^k\right)-A{x}^k\right)}\right)\right\}$ (1.6)

所以问题转化为最小化问题 $\min f\left(x\right)$。我们的目标是找到一个解 $x^{*}$，使在此点的目标函数的梯度满足 $\nabla f\left(x^{*}\right)<\epsilon$。根据上面的替代函数，我们可以推导出非线性分裂可行问题的目标函数和梯度函数分别为

$f\left(x\right)=\frac{1}{2}{\Vert x-P_C\left(x\right)\Vert }^2+\frac{1}{2}{\Vert F\left(x\right)-P_Q\left(F\left(x\right)\right)\Vert }^2$ (1.7)

$\nabla f\left(x\right)=\left(x-P_{C}x\right)+\nabla F{\left(x\right)}^{\text{T}}\left(F\left(x\right)-P_Q\left(F\left(x\right)\right)\right)$ (1.8)

其中 $\nabla F\left(x\right):R^m\to R^n$ 是映射F在x点处的雅可比矩阵。

李志宝 [15] 提出了一种自适应投影类型方法，其文章算法基于投影梯度法加上自适应的技术求解非线性分裂可行问题。作者基于一阶信息使用了新的步长搜索，迭代公式为 $x^{k+1}=P_{\Omega }\left[x^k-{\beta }_{k+1}\nabla p\left(x^k\right)\right]$，其中步长 ${\beta }_{k+1}$ 满足

${\beta }_{k+1}{\Vert \nabla f\left(x^k\right)-\nabla f\left(x^{k+1}\right)\Vert }^2\le \left(2-\delta \right){\left(x^k-x^{k+1}\right)}^{\text{T}}\left(\nabla f\left(x^k\right)-\nabla f\left(x^{k+1}\right)\right)$

2014年A. Gibali [16] 提出了一种求解非线性分裂可行问题的信赖域方法也展现出了良好的收敛性与数值实验结果。

本文，我们给出非线性分裂可行问题的一个条件梯度算法，我们考虑等价的无约束优化问题

$\min f\left(x\right)=\frac{1}{2}{\Vert x-P_C\left(x\right)\Vert }^2+\frac{1}{2}{\Vert F\left(x\right)-P_Q\left(F\left(x\right)\right)\Vert }^2$

该算法将求解迭代方向转化成求解一线性子问题，并以线搜索得到的步长作为凸因子，当前方向与上一步迭代点的凸组合作为新的迭代点。算法在迭代的更新步中不使用投影，并且得到的解有较好的稀疏性和低秩性。我们获得了算法的收敛性并给出数值实验对比分析x了本文的算法与相关算法在同一算例下的表现情况，得到了良好的结果。

2. 预备知识

2.1. 欧式空间上的投影的性质

设向量x为 $R_n$ 空间的任意向量，内积定义为 $\langle ,\rangle$，x的范数定义为 $\Vert x\Vert =\sqrt{\langle x,x\rangle }$。向量x向C集上的投影用 $P_C$ 表示，定义式为

$P_{C}x=\arg \underset{y\in C}{\min }\Vert y-x\Vert ,\forall y\in C$ (2.1)

定理2.1 [17] $P_C$ 与 $I-P_C$ 是非扩张映射，有以下结论成立：

$\Vert P_{C}x-P_{C}y\Vert \le \Vert x-y\Vert$ (2.2)

$\Vert \left(I-P_C\right)\left(x-y\right)\Vert \le \Vert x-y\Vert$ (2.3)

其中 $x,y\in C$。

若算子T满足 $FixT=\left\{x\in D:x=Tx\right\}$，则T是一非扩张算子，称作不动点映射。我们定义

与 $T_Q$ ： $T_C:=\left\{x\in C:x=P_{C}x\right\}$，$T_Q:=\left\{H\left(x\right)=y\in Q:y=P_{Q}y\right\}$ (2.4)

引理2.1 若迭代点 $x^k$ 满足

$x^k\in T_C\cap T_Q$ (2.5)

则说明 $x^k$ 为最优解。

证明：若 $x=x^k\in T_C\cap T_Q$，则 $x^k=P_C{x}^k$ 且 $y^k=F\left(x^k\right)=P_Q\left(x^k\right)$，所以

$f\left(x^k\right)=\frac{1}{2}{\Vert x^k-P_C\left(x^k\right)\Vert }^2+\frac{1}{2}{\Vert F\left(x^k\right)-P_Q\left(F\left(x^k\right)\right)\Vert }^2=0$

显然 $f\left(x\right)\ge 0$，所以当 $x=x^k$ 时，目标函数达到最小值0。根据投影定理，当向量与其投影之间的差为0时，我们可以得出结论：向量 $x^k$ 在集合C的内部或边界上，并且 $F\left(x^k\right)$ 也同时位于集合Q的内部或边界上，即上述两向量满足非线性分裂可行性问题的约束，因此获得了一个解，这意味着向量 $x^k$ 为问题的最优解。

引理2.2 设C是Rn上的闭凸集，则存在

${\left(v-P_C\left(v\right)\right)}^{\text{T}}\left(u-P_C\left(v\right)\right)\le 0$ (2.6)

此不等式关系以二维为例如图1所示如下：

图1. 引理2.2投影不等式示意图

向范数球 $N_D=z\in H:\Vert z-z_0\Vert \le r$ 上投影的解析表达式为

$P_{C}x=\{\begin{array}{l}x,x\in C\\ z_0+\frac{r\left(x-z_0\right)}{\Vert x-z_0\Vert },x\notin C\end{array}$

推论1：假设K是空间X中的闭凸子集。定义任何满足 $x\in X$ 的唯一元素为 ${\pi }_K\left(x\right)\in K$，有

$x-{\pi }_K\left(x\right)={\min }_{y\in K}\Vert x-y\Vert$ (2.7)

其性质为如下变分不等式

$\langle {\pi }_K\left(x\right)-x,{\pi }_K\left(x\right)-y\rangle \le 0,\forall y\in K$ (2.8)

命题1：假设K是一个非空闭凸集，而 $d_K^2$ 是在X上由 $d_K^2\left(x\right)\inf ={\Vert x-y\Vert }^2|y\in K$ 定义的函数。 $d_K^2$ 是连续可微的，其梯度为

$\nabla d_K^2\left(x\right)=2\left(x-{\pi }_K\left(x\right)\right)$ (2.9)

定义等式约束集：

$\mathcal{E}=\left\{x\in R^n:{\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_r\right)}^{\text{T}}x=\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_r\right)\right\}\to \left\{x\in R^n:{\alpha }_i^{\text{T}}x={\beta }_i\right\}$

与不等式约束集：

$\mathcal{J}=\left\{x\in R^n:{\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_r\right)}^{\text{T}}x\le \left({\beta }_1,{\beta }_2,\cdots ,{\beta }_r\right)\right\}\to \left\{x\in R^n:{\alpha }_i^{\text{T}}x\le {\beta }_i\right\}$

向量x向集合 $\mathcal{E}$ 与集合 $\mathcal{J}$ 的投影解析表达式为：

$P_{\mathcal{E}}\left(x\right)=x+\left(b_i-a_i^{\text{T}}x\right)\frac{a_i}{{\Vert a_i\Vert }^2}$

$P_{\mathcal{J}}\left(x\right)=x+\min \left\{0,b_i-a_i^{\text{T}}x\right\}\frac{a_i}{{\Vert a_i\Vert }^2}$

在本文中，我们还将考虑稀疏约束下的非线性分裂可行问题 [18]

$S_r:=\left\{x\in R^t|{\Vert x\Vert }_0\le r,0\le r\le t\right\}$

定义向量x的指标集 $I\left(x\right)=\left\{i_1,i_2,\cdots ,i_r\right\}\subseteq \left\{1,2,\cdots ,n\right\}$，例如 $x=\left\{x_{I_1},x_{I_2},\cdots ,x_{I_n}\right\}$，则其指标集为 $I\left(x\right)=\left\{I_2,I_4,\cdots ,I_{2n}\right\}$，$x=\left\{x_{I_2},x_{I_4},\cdots ,x_{I_{2n}}\right\}$，且有

${\max }_{i\notin I\left(x\right)}\Vert x_i\Vert \le \underset{i\in I\left(x\right)}{\min }\Vert x_i\Vert$

定义x在稀疏集 $S_r$ 上的投影为 $P_{S_r}\left(x\right)$，根据投影的性质，记 $M_s\left(\Vert x\Vert \right)$ 为向量x中最大的分量，我们有：

$P_{C_S}\left(x\right)=\{\begin{array}{l}{x}_i{x}_i>M_s\left(\left|x\right|\right)\\ 0或x_i{x}_i=M_s\left(\left|x\right|\right)\\ 0x_i<M_s\left(\left|x\right|\right)\end{array}$

由于0范数约束问题是一个NP难问题，我们使用范数和范数的最优凸近似 [19] 在问题(1.1)中，如果集合C是一稀疏集合 $C_s=\left\{x\in R^n|{\Vert x\Vert }_0\le r\right\}$，则其可以与下述问题等价

$\underset{x}{\min }f(\; x\; )$

$\text{s}\text{.t}\text{.}{\Vert x\Vert }_1\le \mu$

2.2. 条件梯度法

Marguerite Frank和Philip Wolfe [20] 在1956年首次提出条件梯度法，方法的思想是在每一次迭代中通过求解一个线性的子问题得到最优的下降方向，在更新步时与上一步的迭代点凸组合得到下一步的迭代点。条件梯度法由于其非常好的稀疏性与低秩性，在大数据与机器学习热门的近几年里收货了极大的关注 [21] [22]，其算法描述如下：

${\max }_{\lambda }h(\; \lambda \; )$

$\text{s}\text{.t}\text{.}\lambda \in Q$

在当前迭代点 $x^k$ 处，该算法考虑目标函数的线性化，并求解出该线性函数的最小值(定义域相同)。通过LO oracle解决线性子问题，也可根据不同的问题适当调整子问题的处理方法，从而大大提高了整个算法的灵活性。

以下两个图是条件梯度法的迭代过程的示意图。图2左边表示集合C中的迭代过程。每个迭代方向都来自当前迭代向量，沿着C集合的边界向量的方向。步长可以控制每个迭代向量与边界向量之间的距离以实现更合适的下降。图2右侧图是在分裂可行问题下的，我们的解集是图中的集C和集Q的交集。使用FW方法解决分裂可行问题是非常行之有效的。

图2. 条件梯度法迭代示意图

2.3. 对偶间隙

定义目标函数的对偶间隙 [23] [24]

$\mathcal{G}\left(x\right)=\max \langle x-s,\nabla f\left(x\right)\rangle ,s\in C$ (2.10)

对偶间隙是作为算法每次迭代的副产品，无需进行额外的计算。同时我们记在 $x^k$ 处的下降方向为 $d_k^{FW}=s^k-x^k$。

定义2.1 (下降方向)若目标函数在集合C上是连续可微的，则 $d_k^{FW}$ 为下降方向的充要条件是 $G\left(x^k\right)<0$，并且当 $G\left(x^k\right)=0$ 时，目标函数达到最小值， $x^k$ 为最优解。

证明：由于

$s^{\left(k\right)}:=\arg \max \langle s,\nabla f\left(x^{\left(k\right)}\right)\rangle$

有 $s^k$ 与 $x^k$ 处梯度内积最小，所以

$\langle s^k,\nabla f\left(x^k\right)\rangle \le \langle x^k,\nabla f\left(x^k\right)\rangle$

$\langle s^k,\nabla f\left(x^k\right)\rangle -\langle x^k,\nabla f\left(x^k\right)\rangle \le 0$

$\langle s^k-x^k,\nabla f\left(x^k\right)\rangle \le 0$

$\nabla f{\left(x^k\right)}^{\text{T}}{d}_k^{FW}\le 0$

因此，可以看出， $d_k^{FW}$ 是下降方向，当且仅当 $x^k$ 是最优解时，等号成立。

引理2.4 为一收敛序列，并且有 ${\lim }_{i\to \infty }\mathcal{G}\left(x^i\right)=0$。

证明略。

3. 算法

以下是条件梯度法解决非线性分裂可行问题的算法步骤：

引理3.1 对于 $\forall x,x^i\in C$，$\langle \nabla f\left(x^i\right),x^{*}\rangle \le \langle \nabla f\left(x^i\right),x\rangle$。

证明：目标函数 $f\left(x\right)$ 在 $x^{*}$ 处与任意 $x^r\in C$ 的泰勒展开式为

$f\left(x^{*}\right)=f\left(x^i\right)+\nabla f\left(x^i\right)\left(x^{*}-x^i\right)+\mathcal{O}\left(x\right)=0$

与

$f\left(x^r\right)=f\left(x^i\right)+\nabla f\left(x^i\right)\left(x^r-x^i\right)+\mathcal{O}\left(x\right)={ϵ}_r$

其中 ${ϵ}_r>0$，因此我们有

$\nabla f\left(x^i\right)\left(x^{*}-x^i\right)<\nabla f\left(x^i\right)\left(x^r-x^i\right)$

所以

$\langle \nabla f\left(x^i\right),x^{*}\rangle \le \langle \nabla f\left(x^i\right),x^r\rangle ,\forall x^r,x^i\in C$

证毕。我们自然推导出以下引理

引理3.2 $f\left(x^k\right)-f\left(x^{k-1}\right)<0$

证明

$\begin{array}{c}f\left(x^k\right)-f\left(x^{k-1}\right)=f\left(x^i\right)+\nabla f\left(x^i\right)\left(x^k-x^i\right)-\left[f\left(x^i\right)+\nabla f\left(x^i\right)\left(x^{k-1}-x^i\right)\right]\\ =\nabla f\left(x^i\right)\left(x^k-x^i\right)-\nabla f\left(x^i\right)\left(x^{k-1}-x^i\right)\\ =\langle \nabla f\left(x^i\right),x^k\rangle -\langle \nabla f\left(x^i\right),x^{k-1}\rangle \\ =\nabla f\left(x^i\right)\left(x^k-x^{k-1}\right)={\lambda }_{k-1}\nabla f\left(x^i\right)d_{k-1}^{FW}\le 0\end{array}$

其中 ${\lambda }_{k-1}\in \left(0,1\right)$。当且仅当 $x^i$ 是最优解时，等号成立。

4. 数值实验

4.1. 数值结果

论文应采在本节中我们将通过如下的例子对算法2进行数值实验，并与其他的算法进行对比。

考虑非线性分裂可行问题：

$F\left(x\right):=\left[x_1,x_2,2x_1^2+2x_2^2+1\right]$

非空闭凸集 $C_{NL},C_S$ 与 $Q_{NL}$ 分别是：

(问题1) $C_{NL}:=\left\{x\in R^2|{\left|x\right|}_2\le 1\right\},Q_{NL}:=\left\{y\in R^3|{\left|y-d\right|}_2\le 1\right\}$

(问题2) $C_S:=\left\{x\in R^2|{\left|x\right|}_0\le 1\right\},Q_{NL}:=\left\{y\in R^3|{\left|y-d\right|}_2\le 1\right\}$

其中 $d={\left[0,1.8,3\right]}^{\text{T}}$。

本文的数值实验是使用MATLAB2018a在同一台计算机上进行的。图3是针对上述问题，三种算法的迭代路径示意图。对于问题1，本文比较了最大值最小值法(MM)，传统投影梯度算法(Projected Gradient Method)和本文提出的算法2在不同初始点的迭代次数与计算时间的差异。从表1可以看出，与其他两种算法相比，MM法需要更多的迭代才能到达求解区域。与算法2相比，PGM的迭代路径明显更“曲折”。图4中迭代次数与目标函数值的比较图。可以看出，算法2，PGM和MM方法可以使目标函数值在先前的迭代中迅速减小到较小的值，但是MM需要更多的迭代次数，PGM在迭代的上一步中，函数值将反弹，算法2迅速收敛到0，从而获得良好的结果。表1列出了一些特殊的初

图3. 不同初始点下三种算法的迭代路径比较图

始点对三种算法的迭代次数，迭代时间和最终目标函数值的影响。可以看出，不同的初始点对MM方法有很大的影响，PGM也显示出一定的波动性。在这种情况下，本文提出的算法2具有良好的性能，影响最小。

表1. 较特殊初始点的三种方法求解效率比较表

图4. 不同算法在迭代过程中目标函数值的变化情况

4.2. 不同步长选取策略

算法中步长的选取方式有所不同 [25]，本文主要考虑如下求解步长的方法

1) 经典条件梯度法步长： ${\lambda }_{\omega }=2/k+2$

2) 简单平均步长： ${\lambda }_{\alpha }=1/k+1$

3) 常数步长： ${\lambda }_c=\hat{\lambda }$

4) 线搜索步长(算法2)： ${\lambda }_{LS}$

5) 最小步长选取法： ${\lambda }_{\min }=\left\{1,\frac{{\left(x^k-s^k\right)}^{\text{T}}\nabla f\left(x^k\right)}{2{\Vert x^k-s^k\Vert }^2}\right\}$

我们分析了五种步长选择策略在非线性分裂可行问题中的表现情况，不同步长下的迭代次数比较表如表2所示，图5则对比了迭代过程中目标函数的变化过程情况。

表2. 不同步长下算法各项数值的比较表

图5. 不同步长下迭代过程示意图

5. 结论

本文主要研究利用带线搜索的条件梯度法求解非线性分裂可行问题，根据理论和数值实验的结果，算法在迭代次数、运算时间和实现的可行性方面，优于其他投影类算法，并且灵活性更强。本文仅使用MATLAB中的优化工具箱，没有针对本算例使用特定的方法求解，在求解迭代方向时，可以运用不同的方法，以实现对不同问题的更有效解决方案。对于步长的线搜索，本文使用Armijo线搜索。在以后的研究中，对于更复杂的高纬度问题，可尝试使用非单调线搜索来避免迭代向量非常接近集合的交点时出现的折线步太多的情况。

文章引用

宇振盛,王子伦. 条件梯度法求解非线性分裂可行问题
Solving Nonlinear Split Feasibility Problem via Conditional Gradient Method[J]. 应用数学进展, 2020, 09(09): 1652-1663. https://doi.org/10.12677/AAM.2020.99192

参考文献