分裂可行性问题 (Split Feasibility Problem, SFP), 最早于 1994 年由 Censor 和 Elfving^[1]提出并命名, 其数学模型是

其中非零矩阵 $A\in \mathbb R^{m\times n},\ C\subset \mathbb R^n$ 和 $Q\subset \mathbb R^m$ 是非空闭凸集. 此类问题在医学放射性治疗、信号处理、图像重建、电力系统等领域中有广泛应^[2⇓-4]. 本文始终假设问题 (1.1) 的解集 $\Gamma$ 非空.

由文献 [3,命题 3.1] 知, 在解集 $\Gamma$ 非空的前提下, SFP 可等价转化成如下约束优化问题

其中, $P_{Q}$ 是 $Q$ 上的正交投影 (即 $P_{Q}(y)=\arg\min\{\frac{1}{2}\|z-y\|^2\ |\ z\in Q\},\ \forall\ y\in \mathbb R^m$), $I$ 是单位矩阵, $\|\cdot\|$ 是 $l_2$ 范数. 根据文献 [5], 此目标函数 $f$ 的梯度为

Byrne^[6]最早提出投影型数值解法 (CQ 算法), 其迭代方案是

其中, $P_{C}$ 是 $C$ 上的正交投影, $\gamma>0$ 为步长. 在 CQ 算法框架下, 步长 $\gamma$ 有三种选取方法: 一是固定步长 $\gamma \in(0,\frac{2}{\|A\|^2})$, 步长 $\gamma$ 的选取依赖于 $\|A\|^2$ (${\|A\|^2}$ 表示矩阵 ${A^{{\top}}}A$ 的最大特征值), 其不容易计算, 其数值过大会影响算法的收敛. 二是动态步长, Yang^[7]提出如下动态步长策略

其中, ${\rho_k}>0$, $\sum\limits_{k=1}^{\infty}{\rho_k}={\infty}$, $\sum\limits_{k=1}^{\infty}{{\rho_k}^2}<{\infty}$. 此算法的全局收敛性要求 $Q$ 是有界集、$A$ 是列满秩矩阵. 为了减弱上述收敛性条件, Lopéz 等^[2]提出了新的动态步长策略

三是基于 Armijo 线搜索产生自适应的步长 $\gamma_k$^[8]

其中, $\mu\in{(0,1)}$, $y_{k}=P_C(x_k-{\gamma_k\nabla f(x_k)})$. 为了得到 $\gamma_k$, 通常需多次计算到 $C$ 和 $Q$ 上的投影以及矩阵向量积, 这显然是耗时的.

CQ 算法是求解 SFP 的经典算法. 此外, 文献中也有很多其他方法可求解 SFP. 例如, Wang 和 Xu^[9]提出了求解 SFP 的循环算法; Qin 和 Wang^[10]提出了 Moudaf 型粘性正则化迭代方法; Dong 等^[11]提出了一个线性化分裂算法; 文献 [12] 介绍了几种求解 SFP 的 Dougla-Rachford 分裂算法; 刘等^[13]提出外推加速线性交替方向乘子法求解 SFP.

根据文献 [3], 在 $\Gamma\neq \emptyset$ 的前提下, SFP 的解可等价转化为凸约束非线性单调方程组

的解. 本文致力于从凸约束非线性单调方程组 (1.3) 的视角寻找求解和设计 SFP 的新型算法. 这一思想在目前的文献中未见报道.

目前, 求解一般的约束非线性单调方程组的方法有很多, 如谱梯度投影法^[14,15]、共轭梯度投影法^[16,17]、谱共轭梯度投影法^[18,19]、三项共轭梯度投影法^[20⇓-22]等.

1964 年, Polyak^[23]首次提出了惯性外推技术. 大量的数值实验表明, 惯性外推技术可加快收敛速度. 例如, 针对 SFP, Sahu^[24]和 Suantai 等^[25]在算法中加入惯性外推技术, 有效的提高了算法的性能, 并具有很好的收敛性. 2022 年, Jian 等^[26]提出了非线性伪单调方程组的一簇惯性无导数投影法. 最近, 基于凸组合技术, Jiang 等^[27]提出了求解无约束优化的一簇混合共轭梯度法. 基于文献 [26] 的惯性加速技术和文献 [27] 的凸组合及混合技术, 本文提出一个新型惯性共轭梯度投影算法 (ICGPM). 此算法不用计算矩阵 ${A^\top}A$ 的最大特征值和复杂的多次投影 $P_C$. 在合适假设下, 所提方法具有全局收敛性和 Q-线性收敛率. 数值实验验证了所提方法的有效性.

本文余下部分结构如下: 第 2 节给出算法的具体步骤; 第 3 节证明其全局收敛性和 $Q$-线性收敛率; 第 4 节进行数值实验, 验证所提方法的有效性; 第 5 节给出结论及展望.

记 $\nabla f_k:=\nabla f(x_k)$. 下面给出本文的算法.

算法 1 ICGPM.

步骤 0 任取初始点 $x_0\in C$, 选取参数 $\sigma,\varsigma,\varepsilon>0,\ 0\leq \alpha<1,\ 0<\rho<1,\ 0<\gamma<2,\ \mu>1,$ 置 $x_{-1}:=x_0,\ k:=0$.

步骤 1 若 $\|\nabla f_k\|\leq \varepsilon$, 则停止; 否则, 计算惯性外推步长

并计算惯性步

若 $\|\nabla f(y_k)\|\leq\varepsilon$, 则终止; 否则执行步骤2.

步骤 2 计算搜索方向

其中

凸组合参数 $\theta_k\in[0,1]$. 若 $\|d_k\|\leq\varepsilon$, 则终止; 否则执行步骤3.

步骤 3 采用 Armijo 型线搜索求步长, 即 $t_k=\varsigma{\rho}^{i_k}$, 其中 $i_k$ 是满足下式的最小非负整数 $i$

令 $z_k=y_k+t_kd_k.$

步骤 4 计算下一个迭代点 $x_{k+1}$

其中

置 $k:=k+1$, 转步骤1.

注 2.1 (1) 根据算法 1 易知 $\{x_{k}\}\subset C$, 且在每次迭代中, 算法仅需计算一次到 $C$ 上的投影. (2.3) 式可等价表述为“凸组合”

(2) 在 (2.5) 式中, 参数 $\gamma$ 被称为松弛因子, 其恰当的取值可加速投影收缩方法的数值收敛.

注 2.2 与文献 [26] 相比, 算法 1 采用了不同的搜索方向以及不同的线搜索. 理论上, 除了全局收敛性以外, 本文进一步分析了算法 1 的局部收敛率. 与文献 [27] 相比, 本文设计的搜索方向不需要重启, 且自动满足充分下降性质.

注 2.3 由 (2.1) 式知 $\alpha_k\in[\alpha],$ 且 $\{\alpha_k\}$ 满足

注 2.4 Armijo 型线搜索 (2.4) 式取自文献 [14], 其目的是构造分离超平面. 事实上, 由 (2.4) 式知

进而, 由引理 3.2 得

根据 $\nabla f$ 的单调性 (见命题 3.1(ii)) 知

而 $\nabla f(x^*)=0$ (见命题 3.2), 于是有

由 (2.8) 式和 (2.9) 式知, 超平面

严格分离 $y_k$ 与问题 (1.1) 的解集 $\Gamma$.

本节分析算法 1 的全局收敛性和线性收敛率. 为此, 给出以下命题及引理.

命题 3.1 设函数 $f:\mathbb R^n\rightarrow [0,+\infty)$ 由 (1.2) 式定义, 则

(i) ^[5]$f$ 是连续可微的凸函数, 且其梯度为

(ii) ^[5] 映射 $\nabla f$ 在 $\mathbb R^n$ 上单调, 即 $\nabla f$ 满足

(iii) [4,引理 8.1] 映射 $\nabla f$ 在 $\mathbb R^n$ 上是 $L$-Lipschitz 连续的, 且常数 $L=\|A\|^2>0$, 即

命题 3.2 [3,命题 3.1] 设 SFP 的解集 $\Gamma\neq\emptyset$, 则以下结论等价

(i) $x^*\in\Gamma$;

(ii) $x^*\in C$ 且 $f(x^*)=0$;

(iii) $x^*\in C$ 且 $\nabla f(x^*)=0$.

命题 3.3 对任意 $a,b,c,d\in \mathbb R^{n}$, 有

命题 3.4^[28] 设 $\Omega$ 是 $\mathbb{R}^{n}$ 中非空闭凸子集, 则对任意 $x\in \mathbb{R}^{n}$ 及 $z\in\Omega$, 有

引理 3.1 设 $\{d_k\}$ 是由算法 1 产生的序列, 则对任意 $k\geq1$, 有

进而, 对任意 $k\geq0$, 搜索方向 $d_k$ 满足充分下降条件

及信赖域性质

证 由 (2.3) 式知

其中, $\omega_k\in[{\pi}]$ 是 $\nabla f(y_{k-1})$ 与 $\nabla f(y_{k})$ 的夹角. 由此结合 $\theta_k \in [0,1]$ 知 $\beta_k \geq0$. 进一步,

故 (3.3) 式得证.

下面证明 (3.4) 式. 当 $k=0$ 时, $d_0=-\nabla f(y_0)$, 此时, (3.4) 式显然成立.

当 $k\geq1$ 时, 在 $d_k=-\nabla f(y_k)+\beta_kd_{k-1}$ 两边同时左乘 $\nabla f(y_k)^{\top}$ 得

其中, 上面第一个不等式利用了柯西-施瓦茨不等式及 $\beta_k\geq 0$, 第二个不等式利用了 (3.3) 式. 因此, (3.4) 式得证.

接下来, 证明 (3.5) 式. 当 $k=0$ 时, $d_0=-\nabla f(y_0)$, 则 (3.3) 式成立.

当 $k\geq1$ 时, $d_k=-\nabla f(y_k)+\beta_kd_{k-1}$. 于是, 利用三角不等式及 (3.3) 式得

故 (3.5) 式得证. 证毕.

引理 3.2 算法 1 步骤 3 的 Armijo 型线搜索可在有限次计算后完成, 从而算法 1 是适定的.

证 反证法. 假设存在某个非负整数 $k_0$, 对任意非负整数 $i$, (2.4) 式均不成立, 即

令 $i\rightarrow\infty$, 再结合 $\nabla f$ 的连续性, 可知 $\nabla f(y_{k_0})^\top d_{k_0}\geq0$. 这与 (3.4) 式矛盾, 故引理成立.

引理 3.3 设 SFP 的解集 $\Gamma\neq\emptyset$, 则对任意 $x^*\in\Gamma$, 算法 1 产生的迭代序列 $\{x_k\},\ \{y_k\}$ 和 $\{z_k\}$ 满足

证 由 $x^*\in\Gamma$ 及命题 3.2 得

此结合 (3.1) 式知

根据 $z_k$ 定义及 (2.8) 式得

于是, 由 (2.5) 式和命题 3.4 知

其中, 上面第三个不等式利用了 (2.6) 式, (3.7) 式及 (3.8) 式, 最后一个不等式由 $\xi_k$ 定义和 (3.8) 式得到. 由 (2.2) 式知

其中, 上面第三个等式由命题 3.3 得到, 最后一个不等式利用了 $\alpha_k\in[\alpha]\subset[0,1)$. 结合 (3.9) 式及 (3.10) 式可得 (3.6) 式. 证毕.

引理 3.4^{([29], [30,引理 7])} 设非负数列 $\{\chi_k\}$, $\{\eta_k\}$, $\{\alpha_k\}$ 和 $\{\upsilon_k\}$ 满足 $\chi_{0}=\chi_{-1}$, $\alpha_k\in[\alpha]\subset[0,1)$, 且

若 $\sum\limits_{k=0}^{\infty}\eta_k<\infty$, 则

(i) $\lim\limits_{k \rightarrow \infty}\chi_{k}$ 存在;

(ii) $\sum\limits_{k=0}^{\infty}\upsilon_k<\infty$.

定理 3.1 设 $\{x_k\},\ \{y_k\}$, $\{z_k\}$ 和 $\{\alpha_k\}$ 由算法 1 产生, 且 $x^*\in\Gamma$, 则

(i) $\lim\limits_{k \rightarrow \infty}\|x_k-x^*\|$ 存在;

(ii) $\sum\limits_{k=0}^{\infty}\frac{\|z_k-y_k\|^4}{\|\nabla f(z_k)\|^2}<\infty$.

证 记

于是, 根据 (3.6) 式得

由 $\gamma\in (0,2)$, $\sigma>0$, (2.7) 式及引理 3.4 知定理成立.

推论 3.1 若定理 3.1 的条件成立, 则序列 $\{x_k\}$ 和 $\{y_k\}$ 均有界, 且

证 由 (2.7) 式得

此结合 (2.2) 式和 $\alpha_k\in[\alpha]\subset[0,1)$, 推得

另外, 根据定理 3.1(i) 知序列 $\{x_k\}$ 有界. 再结合 (3.12) 式知序列 $\{y_k\}$ 也有界. 最后, 由定理 3.1}(ii) 立得 (3.11) 式成立. 证毕.

下面给出本文的全局收敛性定理

定理 3.2 由算法 1 产生的迭代序列 $\{x_k\}$, $\{y_k\}$ 和 $\{z_k\}$ 均收敛到 SFP (1.1) 的解, 且 $\lim\limits_{k\rightarrow \infty}x_{k}=\lim\limits_{k\rightarrow \infty}y_{k}=\lim\limits_{k\rightarrow \infty}z_{k}$.

证 首先, 证 $\lim\limits_{k\rightarrow\infty}\|z_k-y_k\|=\lim\limits_{k\rightarrow\infty}t_k\|d_k\|=0$. 依据推论 3.1 知 $\{x_k\}$ 和 $\{y_k\}$ 有界. 此结合 (3.5) 式及 $\nabla f$ 连续性知 $\{d_k\}$ 有界. 于是, 由 $t_k\in(0,\varsigma]$ 及 $z_k$ 定义易知 $\{z_k\}$ 有界. 因此, 由 $\nabla f$ 连续性可知, 存在常数 $M_0>0$, 使得

此结合 (3.11) 式, 有

因此, 由 $z_k$ 定义知

其次, 证明 $\underset{k\rightarrow\infty}{\lim}\|\nabla f(y_k)\|=0$. 反证法. 假设 $\underset{k\rightarrow\infty}{\limsup}\|\nabla f(y_k)\|\neq0$, 则存在常数 $\varepsilon_{0}>0$ 及无穷子列 $\mathcal{K}$, 使得

由 (3.4) 式知

进而, 由 $\underset{k\rightarrow\infty}{\lim}t_k\|d_k\|=0$ 得 $\lim\limits_{k\in \mathcal{K}}t_k=0$. 根据 $\{d_k\}$ 和 $\{y_k\}$ 的有界性, 不妨设

在 (3.4) 式中令 $k\in \mathcal{K}$, 取极限, 并利用 (3.15) 式, (3.17) 式和 $\nabla f$ 的连续性, 得

由算法 1步骤 3 的 Armijo 型线搜索准则知, 对充分大的 $k\in \mathcal{K}$, 有

同理, 对上述不等式取极限, 并利用 $\lim\limits_{k\in \mathcal{K}}t_k=0$, (3.17) 式与 $\nabla f$ 的连续性, 得 $-\nabla f(\bar{y})^\top \bar{d}\leq0$. 此与 (3.18) 式矛盾. 故 $\underset{k\rightarrow\infty}{\limsup}\|\nabla f(y_k)\|=0$, 进而

成立.

最后, 验证 $\{x_k\}$, $\{y_k\}$, $\{z_k\}$ 均收敛到 SFP (1.1) 的解. 取有界序列 $\{y_k\}$ 的一个收敛子列 $y_{k_i}\rightarrow\bar{y},\ i\rightarrow\infty$. 于是, 根据 (3.13) 式, (3.19) 式和 $\nabla f$ 的连续性, 有

另一方面, 由(3.12) 式易得 $x_{k_i}\rightarrow\bar{y}\ (i\rightarrow\infty)$. 因此, 由 $\{x_k\}\subset C$ 及 $C$ 的闭性知, $\bar{y}\in C$. 此结合 (3.20) 式的第二个关系式及 (1.3) 式知 $\bar{y}\in\Gamma$. 在定理 3.1 中, 令 $x^*:=\bar{y}$, 利用 $\lim\limits_{i\rightarrow\infty}x_{k_i}=\bar{y}$, 得

此进一步表明整个序列 $\{x_k\}$ 收敛到 $\bar{y}\in\Gamma$. 于是, 由 (3.12) 式和 (3.14) 式推出序列 $\{y_k\}$ 和 $\{z_k\}$ 也收敛到 $\bar{y}\in\Gamma$. 证毕.

基于定理 3.2, 在如下假设条件下, 我们进一步讨论算法 1 的线性收敛速度.

假设 3.1 对序列 $\{y_k\}$ 的极限 $\bar{y}\in \Gamma$, $\nabla f$ 在 $\bar{y}$ 附近满足局部误差界条件, 即存在两个常数 $l>0,\ \epsilon>0$, 使得对任意 $y\in N(\bar{y},\epsilon)$, 有

其中 ${\rm dist}(y,\Gamma)$ 表示点 $y$ 到解集 $\Gamma$ 的距离, 邻域 $N(\bar{y},\epsilon):=\{y\in \mathbb R^n\ |\ \|y-\bar{y}\|<\epsilon\}$.

注 3.1 若 $\nabla f$ 是 $l$-强单调的, 即

则假设 3.1 成立. 事实上, 由柯西-施瓦茨不等式知

于是有

令 $x:=\bar{y}\in \Gamma$, 则

下面引理给出步长序列 $\{t_k\}$ 的一个正的下界.

引理 3.5 对于算法 1 产生的步长序列 $\{t_k\}$, 有

证 令 $\overline{t}_k=\rho^{-1}t_k$, 对于 $t_{k}<\varsigma$, 根据 Armijo 型线搜索准则知,

令 $\overline{z}_k=y_k+\overline{t}_k d_k$, 结合 (3.2) 式, (3.4) 式和 (3.5) 式, 得

故 $t_k\geq\frac{\rho(1-\frac{1}{\mu})}{(L+\sigma)(1+\frac{1}{\mu})^2}$. 证毕.

借助于文献 [31] 关于线性收敛率的分析技巧, 下面建立算法 1 的线性收敛率.

定理 3.3 设假设 3.1 成立, 序列 $\{x_k\}$ 由算法 1 产生, 且记 $a=\underset{k\rightarrow\infty}{\lim\sup}\frac{\frac{1}{k}}{{\rm dist}(x_k,\Gamma)}$.

(1) 若 $0<a<\infty$, 则至少存在一个无限指标集 $\mathcal{K}$, 使得 ${\rm dist}(x_k,\Gamma)=O(\frac{1}{k}),\ k\in \mathcal{K}$;

(2) 若 $a=\infty$, 则至少存在一个无限指标集 $\mathcal{K}$, 使得 ${\rm dist}(x_k,\Gamma)=o(\frac{1}{k}),\ k\in \mathcal{K}$;

(3) 若 $a=0$, 则序列 $\{{\rm dist}(x_k,\Gamma)\}$ Q-线性收敛于 0, 即 $\underset{k\rightarrow\infty}{\lim\sup}\frac{{\rm dist}(x_{k+1},\Gamma)}{{\rm dist}(x_k,\Gamma)}<1.$

证 结论 (1) 和结论 (2) 显然成立, 现证明结论 (3).

设 $h_k$ 是 $\Gamma$ 中最接近 $y_k$ 的点, 即 $\|y_k-h_k\|={\rm dist}(y_k,\Gamma)$, 则由 (3.9) 式知

另一方面, 根据 $\nabla f$ 的 Lipschitz 连续性, (3.5) 式, $t_k\in(0,\varsigma]$, 并注意到 $\nabla f(h_k)=0$, 得

根据 $z_k$ 定义, (3.16) 式, (3.21) 式及 (3.22) 式知, 当 $k$ 充分大时, 有

结合 (3.23) 式, (3.24) 式和 (3.25) 式, 当 $k$ 充分大时, 可得

其中

由 (2.1) 式知 $\alpha_k\|x_k-x_{k-1}\|^2\leq\frac{1}{k^2}$. 此结合 $\alpha_k\in[\alpha]\subset[0,1)$ 得

进而,

其中, $h_{x_k}$ 是 $\Gamma$ 中最接近 $x_k$ 的点. 因此

此表明若 $\underset{k\rightarrow\infty}{\lim\sup}\frac{\frac{1}{k}}{{\rm dist}(x_k,\Gamma)}=0$, 则结论 (3) 成立. 证毕.

本节将通过 5 个例子检验算法 1 (ICGPM) 的数值效果. 所有程序代码利用 MATLAB R2017a 编写, 软件运行环境为 Windows 10 家庭中文版, Intel(R) Core(TM) i5-6200U, CPU 2.40GHz, 4GB RAM. 所有程序的终止准则为

本文将所提算法 ICGPM 和 IPCM^[25]相比较. 测试问题如下

问题 1 求 $x\in C:=\{x\in\mathbb R^3\ |\ \|x\|\leq2\}$, 使得 $Ax\in Q$, 其中

问题 2 求 $x\in C:=\{x\in \mathbb R^3\ |\ x_1+x_2+x_3\leq3\}$, 使得 $Ax\in Q$, 其中

问题 3 求 $x\in C:=\{x\in\mathbb R^4\ |\ \|x\|\leq2\}$, 使得 $Ax\in Q$, 其中

问题 4 求 $x\in C:=\{x\in\mathbb R^4\ |\ \|x\|\leq2\}$, 使得 $Ax\in Q$, 其中

问题 5 求 $x\in C:=\{x\in \mathbb R^3\ |\ x_1+x_2+x_3\leq1\}$, 使得 $Ax\in Q$, 其中

注 4.1 对于上述问题, 算法 1 步骤 4 中投影的计算参考文献 [32,例 29.20] 和文献 [33,命题 2.2.1].

经过反复调试参数, 本文选取两个算法较好的参数取值进行数值比较. 下面给出两个算法的具体参数设置

ICGPM: $\varsigma= 1.15, \sigma = 0.0005, \rho= 0.4, \gamma= 1.69, \theta_k\equiv 0.5, \mu=3, \alpha=0.63.$

IPCM: $\gamma= 1.2, \sigma= 0.015$,$\mu= 0.4, \rho= 0.2, \alpha=0.8.$

在数值实验中, 每个测试问题所采用的初始点 (IP) 如下

问题 1: ${x}_0^1=(0,0,0)^\top, {x}_0^2=(1,1,1)^\top, {x}_0^3=(0,-2,0)^\top, {x}_0^4=(-1,1,-1)^\top.$

问题 2: ${x}_0^5=(2,-4,3)^\top, {x}_0^6=(1,1,1)^\top, {x}_0^7=(10,8,2)^\top, {x}_0^8=(-1,-2,3)^\top.$

问题 3: ${x}_0^9=(0,0,0,0)^\top, {x}_0^{10}=(6,4,20,6)^\top, {x}_0^{11}=(1,5,6,-2)^\top, {x}_0^{12}=(5,-1,10,8)^\top.$

问题 4: ${x}_0^{13}=(-30,-20,40,-5)^\top, {x}_0^{14}=(0,-3,-10,-5)^\top, {x}_0^{15}=(-10,0,10,2)^\top, {x}_0^{16}=(5,20,28,35)^\top.$

问题 5: ${x}_0^{17}=(-2,-4,3)^\top, {x}_0^{18}=(-10,8,-7)^\top, {x}_0^{19}=(-8,1,0)^\top, {x}_0^{20}=(9,5,-20)^\top.$

在表1中, “Itr”表示迭代次数, “NF”表示 $\nabla f$ 的计算次数, $\|\nabla f^\ast\|$ 表示 $\|\nabla f\|$ 在程序终止时的值. 根据表1, 对于所测试的例子, 可得如下结论

(1) 本文所提算法 1 可有效求解 SFP;

(2) 基于问题 1 和问题 4 的数值结果可见, 本文所提算法 ICGPM 明显优于 IPCM;

(3) 从计算时间及解的质量看, 算法 1 迭代次数少、运行时间短, 且能获得更高质量的近似解.

本文基于非线性单调方程组的求解方法, 提出了求解分裂可行性问题的惯性共轭梯度投影法. 在较弱的条件下, 证明了算法的全局收敛性及$Q$-线性收敛率. 数值结果初步验证了算法的有效性和鲁棒性.

注意到, 本文所提算法的主要计算量为投影算子 $P_Q(\cdot)$ 和 $P_C(\cdot)$ 的计算. 当 $Q$ 和 $C$ 较复杂时, 相应投影算子的计算是耗时的和高成本的. 结合本文所提算法思想与松弛技术^[7,34,35], 设计惯性梯度类松弛投影法是一个值得研究的课题.