0 引言

在振动系统、潮流方程等科学与工程计算领域存在许多大规模优化问题^{[1, 2]}，而这些优化问题往往能够转化为非线性方程组问题。因此，研究求解大规模非线性方程组的高效数值算法具有重要的理论价值与实际意义。

本文主要考虑以下非线性方程组问题：

$ F\left( x \right) = 0, x \in {R^n}, $

(1)

其中F:Rⁿ→Rⁿ是连续可微非线性函数。令$ f\left( x \right) = \frac{1}{2}{\left\| {F\left( x \right)} \right\|^2} $，其中‖·‖为Euclidean范数，那么求解非线性方程组问题(1)等价于求解如下无约束优化问题：

$ \min f\left( x \right), x \in {R^n}。$

近年来，求解上述优化问题的常见算法有牛顿法、信赖域法、拟牛顿法、Levenberg-Marquardt算法及其各种变形^[3-8]。在选择合理初始点的前提下，上述算法对于小规模优化问题具有快速收敛和数值效果良好等特点，但在迭代过程中，需要计算和存储相关矩阵信息，给求解大规模优化问题带来一定的局限性。为建立求解大规模优化问题的高效算法体系，研究者提出具有算法简单、计算和存储量低等优点的共轭梯度法^[9-11]。

经典共轭梯度法的一般迭代公式为

$ {x_{k + 1}} = {x_k} + {a_k}{d_k}, k = 0, 1, 2, \ldots , $

其中α_k为由某种线搜索所决定的步长。搜索方向d_k为

$ {d_k} = \left\{ \begin{array}{l} \;\;\;\;\; - {F_k}, \;\;\;\;\;\;\;k = 0\\ - {F_k} + {\beta _k}{d_{k - 1}}, k \ge 1 \end{array} \right., $

其中, β_k为共轭参数，F_k为F(x_k)的简写。

本文基于Abubakar等^[12]提出的修正FR搜索方向，借鉴Yuan等^[13]的凸组合思想，构造凸组合系数如下：

$ {N_k} = \frac{{{{\left\| {{y_{k - 1}}} \right\|}^2}}}{{y_{k - 1}^t\omega _{k - 1}^ * }}, $

其中：$ {y_{k - 1}} = {F_k} - {F_{k - 1}}, \omega _{k - 1}^ * = {\omega _{k - 1}} + \left( {\max \left\{ {0, \frac{{ - \omega _{k - 1}^{\rm{T}}{y_{k - 1}}}}{{{{\left\| {{y_{k - 1}}} \right\|}^2}}}} \right\} + 1} \right){y_{k - 1}}, {\omega _{k - 1}} = {x_k} - {x_{k - 1}} $。

同时，采用Andrei^[14]的加速线搜索技术，提出一个求解大规模非线性方程组问题的加速FR型共轭梯度算法。

1 算法描述与性质

本节主要讨论搜索方向的构建并介绍线搜索技术，同时提出凸组合修正共轭梯度算法。

首先，Abubakar等^[12]在2019年提出一种修正FR共轭梯度法，其搜索方向为

$ \begin{array}{l} {d_k} = \\ \left\{ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - {F_k}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k = 0\\ - {F_k} + \frac{{{{\left\| {{F_k}} \right\|}^2}{\omega _{k - 1}} - F_k^{\rm{T}}{\omega _{k - 1}}{F_k}}}{{\max \left\{ {\mu \left\| {{\omega _{k - 1}}} \right\|\left\| {{F_b}} \right\|, {{\left\| {{F_{b - 1}}} \right\|}^2}} \right\}}}, k \ge 1 \end{array} \right., \end{array} $

其中, ω_k-1 =x_k -x_k-1, μ >0。该搜索方向具备充分下降性和信赖域特征，能有效求解大规模无约束优化问题。基于Yuan等^[13]的凸组合思想，本文构建一个新型的凸组合搜索方向：

$ {d_k} = \left\{ \begin{array}{l} - {F_k}, k = 0\\ - {N_k}{F_k} + \left( {1 - {N_k}} \right)\left( {{{\left\| {{F_k}} \right\|}^2}{\omega _{k - 1}} - F_k^{\rm{T}}{\omega _{k - 1}}{F_k}} \right)/\\ \;\;\;\max \left\{ {2\mu \left\| {{\omega _{k1}}} \right\|\left\| {{F_k}} \right\|, {{\left\| {{F_{k1}}} \right\|}^2}} \right\}, k \ge 1 \end{array} \right., $

(2)

其中，N_k满足如下性质：$ y_{k - 1}^{\rm{T}}\omega _{k - 1}^ * \ge \max \left\{ {y_{k - 1}^{\rm{T}}} \right.{\omega _{k - 1}}, \left. {{{\left\| {{y_{k - 1}}, } \right\|}^2}} \right\} \ge {\left\| {{y_{k - 1}}} \right\|^2} > 0 $，对任意的k都有N_k∈ (0, 1]。

其次，本文通过下述方法计算步长α_k=r^m_k，使得m_k满足下式的最小非负整数，即

$ f({x_k} + {r^m}{d_k}) - f({x_k}){\rm{ }} \le \sigma {({r^m})^2}F_k^{\rm{T}}{d_k}, $

(3)

其中，σ∈(0, 1), r∈(0, 1)。Andrei^[14]研究表明，合理地应用加速线搜索，将有效提高算法的计算效率。于是借鉴于Andrei^[14]的加速线搜索技术思想，对步长α_k做出修正，即

$ \tilde \alpha = {\gamma _k}{\alpha _k}, $

其中$ {\varphi _k} = {\alpha _k}F_k^{\rm{T}}{d_k}, {\theta _k} = - {\alpha _k}y_{k - 1}^{\rm{T}}{d_k}, {\gamma _k} = - \frac{{\varphi k}}{{\theta k}} $。若θ_k>0，则记α_k=ã。

最后，建立求解非线性方程组问题(1)的凸组合加速FR型共轭梯度算法(MMFR)。

步骤1：给定初始点x₀∈Rⁿ，参数ε, σ, r, β, μ∈(0, 1)，令k:=0；

步骤2：若‖ F_k‖≤ε，则算法停止；

步骤3：通过式(2)计算搜索方向d_k；

步骤4：若‖F(x_k+d_k)‖≤β‖F_k‖，则令步长α_k=1，转步骤6，否则转步骤5；

步骤5：通过式(3)决定步长α_k；

步骤6：计算θ_k和ã ，若θ_k>0，则令α_k=ã，否则保持步长α_k不变；

步骤7：(更新步)更新新的迭代点x_k+1=x_k+α_kd_k，令k:=k+1，转步骤2。

为后续证明算法的全局收敛性质，下面分析搜索方向d_k的两个重要性质：充分下降性和信赖域特性。

引理1  算法MMFR产生的序列 d_k 和 F_k 满足以下性质：

$ F_k^{\rm{T}}{d_k} \le - {N_k}{\left\| {{F_k}} \right\|^2}, $

(4)

及

$ \left\| {{d_k}} \right\|{\rm{ }} \le {\tau _k}\left\| {{F_k}} \right\|{\rm{ }}, $

(5)

其中，$ {\tau _k} = {\rm{ }}\left( {{N_k} + (1 - {N_k})\frac{1}{\mu }} \right) $

证明：当k=0时，F₀^Td₀=-‖F₀‖², ‖d₀‖=‖F₀ ‖，显然式(4)(5)成立。当k≥1时，由式(2)左乘F_k^T得

$ \begin{array}{l} F_k^{\rm{T}}{d_k} = - {N_k}{\left\| {{F_k}} \right\|^2} + \\ \left( {1 - {N_k}} \right)\frac{{{{\left\| {{F_k}} \right\|}^2}\omega _{k - 1}^{\rm{T}}{F_k} - F_k^{\rm{T}}{\omega _{k - 1}}F_k^{\rm{T}}{F_k}}}{{{\rm{max }}\left\{ {2\mu \left\| {{\omega _{k - 1}}} \right\|{\rm{ }}\left\| {{F_k}} \right\|{\rm{ }}, {\rm{ }}{{\left\| {{F_{k - 1}}} \right\|}^2}} \right\}}} \le \\ - {\rm{ }}{N_k}{\left\| {{F_k}} \right\|^2}。\end{array} $

此外，由式(2)和Cauchy-Schwartz不等式可知

$ \begin{array}{*{20}{l}} {\left\| {{d_k}} \right\| = \left\| { - {N_k}{F_k} + \left( {1 - {N_k}} \right)} \right. \cdot }\\ \begin{array}{l} \left. {\frac{{{{\left\| {{F_k}} \right\|}^2}{\omega _{k - 1}} - F_k^{\rm{T}}{\omega _{k - 1}}{F_k}}}{{\max \left\{ {2\mu \left\| {{\omega _{k - 1}}} \right\|\left\| {{F_k}} \right\|,{{\left\| {{F_{k - 1}}} \right\|}^2}} \right\}}}} \right\| \le \\ {N_k}\left\| {{F_k}} \right\| + \left( {1 - {N_k}} \right). \end{array}\\ {\frac{{{{\left\| {{F_k}} \right\|}^2}\left\| {{\omega _{k - 1}}} \right\| + \left\| {{F_k}} \right\|\left\| {{\omega _{k - 1}}} \right\|\left\| {{F_k}} \right\|}}{{2\mu \left\| {{\omega _{k - 1}}} \right\|\left\| {{F_k}} \right\|}} \le }\\ {\left( {{N_k} + \left( {1 - {N_k}} \right)\frac{1}{\mu }} \right)\left\| {{F_k}} \right\|,} \end{array} $

综上所述，式(4)和式(5)成立, 引理1得证。

2 全局收敛性分析

为进一步分析算法MMFR的收敛性，本节做如下假设：

假设H

(H1) 函数F(x)在开凸集Ω₁⊆Ω= {x|‖F (x) ‖≤‖F (x₀)‖是连续可微的；

(H2) 函数F (x)的雅可比矩阵为▽F (x)是有界的且为对称正定矩阵，即存在正常数ξ₁≥ξ₂>0，使得有‖▽F (x) ‖≤ξ₁和ξ₂‖p‖²≤p^T▽F (x) p≤ξ₁‖p‖², p∈Rⁿ。

引理2 在假设H条件下，序列{ x_k}、{α_k}、{F_k}和{d_k}是由算法MMFR产生的，那么{‖ F_k‖ }为单调递减的收敛序列且$ \mathop {\lim }\limits_{k \to \infty } \left\| {{a_k}{F_k}} \right\| = 0 $。

证明：由Brown等^[15]的引理3.8可得

f(x_k+r^md_k)-f(x_k)≤σ₁(r^m)²F_k^T▽F_kd_k，利用假设H2的公式，结合上式得f(x_k+r^md_k)-f(x_k)≤σ(r^m)²F_k^Td_k，其中σ=σ₁ξ₁。因此，线搜索是适定的。

由引理1和式(3)可推出

$ f({x_k} + {\alpha _k}{d_k}) - f({x_k}){\rm{ }} \le - \sigma \alpha _k^2{N_k}{\left\| {{F_k}} \right\|^2} \le 0, $

(6)

这说明函数f(x)沿着下降方向d_k是充分下降的。公式(6)结合f(x)的定义可知，对于任意的k，都有‖F_k+1‖≤‖F_k‖。此外，由式(3)和式(4)得出

$ \begin{array}{l} \mathop {\lim }\limits_{k \to \infty } f\left( {{x_k} + {a_k}{d_k}} \right) - f\left( {{x_0}} \right) = \\ \sum\limits_{k = 0}^\infty {f\left( {{x_k} + {a_k}{d_k}} \right) - f\left( {{x_k}} \right) \le - \sigma \sum\limits_{k = 0}^\infty {{N_k}{{\left\| {{a_k}{F_k}} \right\|}^2}} } 。\end{array} $

由假设H1中函数的有界性，再结合上式得

$ \begin{array}{l} \;\;\;\;\;\;\sigma \sum\limits_{k = 0}^\infty {{N_k}{{\left\| {{a_k}{F_k}} \right\|}^2}} \le \mathop {\lim }\limits_{k \to \infty } f\left( {{x_0}} \right) - f({x_k} + \\ {a_k}{d_k}) < \infty , \end{array} $

再结合N_k∈ (0, 1]推得$ \mathop {\lim }\limits_{k \to \infty } f\left( {{x_0}} \right)\left\| {{a_k}{F_k}} \right\| = 0 $ 。

下面给出算法MMFR的全局收敛性定理。

定理1  在假设H条件下，序列{ x_k}、{α_k}、{ F_k}和{d_k}是由算法MMFR产生的，那么有$ \mathop {\lim }\limits_{k \to \infty } \left\| {\nabla f\left( {{x_k}} \right)} \right\| = 0 $ 。

证明:由引理2可知，若$ \mathop {\lim }\limits_{k \to \infty } {a_k} > 0, 则\mathop {\lim }\limits_{k \to \infty } \left\| {{F_k}} \right\| = 0 $，结合▽f(x)=▽F(x)F(x)可知结论成立；若$ \mathop {\lim }\limits_{k \to \infty } {a_k} = 0, 则显然有\mathop {\lim }\limits_{k \to \infty } {a_k} \ne r $，那么令α_k=α_kr^-1，于是有如下不等式成立：

$ f\left( {{x_k} + {{\overline \alpha }_k}{d_k}} \right) - f\left( {{x_k}} \right) > \sigma \overline \alpha _k^2F_k^{\rm{T}}{d_k}。$

(7)

假设结论不成立，即存在正整数ξ，对于任意k，那么有

$ \left\| {\nabla f\left( {{x_k}} \right)} \right\| > \xi 。$

(8)

另外，由假设H1可知，集合为有界集合，则序列{x_k }是有界的，于是可得序列{d_k}也是有界的。不失一般性，设点x^*和d^*分别为序列{ x_k }和{ d_k }的聚点。因此，对式(7)取极限得

$ \nabla f{\left( {{x^ * }} \right)^T}{d^ * } \ge 0。$

同理，对式(6)取极限得

$ \nabla f{\left( {{x^ * }} \right)^{\rm{T}}}{d^ * } \le 0。$

上述两不等式可推出▽f(x^*)^Td^*=0，由于‖d^*‖≠0 (如果‖d^*‖=0，那么由式(4)知‖F(x^*)‖=0)，故‖▽f(x^*)‖=0，这与式(8)矛盾，于是结论成立。

定理1说明序列{ x_k }至少是线性收敛的，下面定理给出算法MMFR具有强收敛性质。

定理2  在假设H条件下，若算法MMFR产生的任一子序列{ x_k }收敛于聚点x^*，则非线性方程组问题(1)的最优解为x^*，进一步有序列{ x_k }整列收敛于x^*。

证明：类似于Yuan^[16]中定理3.4的证明方法，易证结论成立，故省略证明过程。

3 数值试验

本节通过比较算法MMFR、经典FR、三项FR算法在求解大规模非线性方程组问题上的数值结果，验证算法MMFR的有效性与稳定性。

下面给出经典FR算法和三项FR算法的搜索方向，分别为

$ {d_k} = \left\{ \begin{array}{l} \;\;\;\;\;\;\;\;\;\; - {F_k}, \;\;\;\;\;\;\;\;\;k = 0\\ - {F_k} + \frac{{{{\left\| {{F_k}} \right\|}^2}}}{{{{\left\| {{F_{k - 1}}} \right\|}^2}}}{d_{k - 1}}, k \ge 1 \end{array} \right., $

(9)

$ \begin{array}{l} \;\;\;\;\;{d_k}{\rm{ = }}\\ \left\{ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - {F_k}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k = 0\\ - {F_k} + \frac{{{{\left\| {{F_k}} \right\|}^2}}}{{{{\left\| {{F_{k - 1}}} \right\|}^2}}}{\omega _{k - 1}}{\rm{ + }}\frac{{F_k^{\rm{T}}{\omega _{k - 1}}}}{{{{\left\| {{F_{k - 1}}} \right\|}^2}}}{F_k}, k \ge 1 \end{array} \right.。\end{array} $

(10)

在算法MMFR的步骤2中，分别采用式(9)和式(10)产生搜索方向，其余步骤不变，得到的算法记为FR算法和MFR算法。

参数设置：r=0.5, σ=0.068, μ=0.25, β=0.5。程序运行环境：MATLAB(2014a)软件实现，Windows10 (64 bite)，RAM:8 G，CPU 3.60 GHz。

算法终止准则为‖F_k‖≤10^-5或Iter>3000，维数为[4500, 12000, 24000, 30000, 45000]。测试问题的函数名称和初始点^[17]见表 1，数值试验结果如表 2所示，其中Pro (Problem)为问题序号，Dim (Dimension)为维数，Iter (Iterations)为迭代次数，NF(The number of function)为函数F(x)计算次数，Time为程序运行时间(单位：s)。根据迭代次数、函数计算次数和运行时间可以看出，总体上算法MMFR最好，算法MFR其次，算法FR最差(表 2)。

为直观地展示3种算法的性能差异，本节采用性能曲线描绘方法^[18]分别描绘出迭代次数性能图、函数计算次数性能图和运行时间性能图(图 1-3)。由图 1-3可知，算法MMFR总体上比算法MFR和算法FR更优，且具有更好的鲁棒性，因此本文提出的算法MMFR是有效的和鲁棒的。

4 结论

本文在修正FR算法的基础上，结合凸组合思想，构造出一个新的修正搜索方向，并利用加速线搜索技术，提出一个加速FR型共轭梯度算法。新的搜索方向不依赖线搜索, 具有充分下降和信赖域特性，还具有良好的理论性质与数值效果。因计算简单，存储量小，十分适合求解大规模非线性方程组问题。同时，也可尝试将新算法进一步推广到信号恢复等实际应用中。