求解对称非负逆特征值问题的一阶黎曼算法

Riemannian Algorithms of Order One for Solving Symmetric Nonnegative Inverse Eigenvalue Problems

摘要: 基于对称非负逆特征值问题的黎曼优化模型，提出了求解该问题的黎曼梯度下降算法和黎曼共轭梯度算法，并分析了这两种算法的收敛性。通过数值实验，比较了两种算法的收敛效率。数值结果表明，求解高阶的非负逆特征值问题时，黎曼共轭梯度算法的收敛效率要高于黎曼梯度下降算法。

Abstract: The Riemannian gradient descent algorithm and the Riemannian conjugate gradient algorithm are proposed for solving the symmetric nonnegative inverse eigenvalue problems, respectively, to-gether with the convergence analysis. At last, some numerical experiments are provided to illus-trate the convergence of the algorithms, and they show that the Riemannian conjugate gradient al-gorithm is more effective than the Riemannian gradient descent algorithm.

1. 引言

矩阵中所有元素都是非负实数的矩阵称为非负矩阵。非负矩阵理论一直是矩阵理论中重要的研究领域之一，并且在数学和自然科学的其它分支以及社会科学中都被广泛涉及，例如博弈论、Markov链、概率论、数值分析、离散分布、群论、小振荡弹性系统和经济学等。经典的非负逆特征值问题是已知矩阵的全体特征值，求非负矩阵(简称NIEP)。近来，白正简老师等人研究了已知矩阵部分特征对的非负逆特征值问题 [1] [2] [3]。特别地，若所求非负矩阵是对称的，则称为对称非负逆特征值问题(简称SNIEP)；若所求非负矩阵是随机矩阵(行和为一的非负矩阵)，则称为对称随机逆特征值问题(简称StNIEP)。NIEP在物理学、控制设计、结构动力学、博弈论、马尔可夫链、概率算法、数值分析、离散分布，分类数据，群论和经济学等领域有广泛的应用 [4] [5] [6]。自从1937年数学家Kolmogorov提出矩阵非负逆特征值问题以来，一直备受关注。关于NIEP的研究主要是两个方面，一方面关于解的理论研究，包括解的存在性和唯一性，另一方面算法研究，即研究求解各类非负逆特征值问题的数值算法研究。关于第一个方面的研究，虽然取得了一些重要的研究结果，至今仍未彻底解决。1949年数学家Sulemanova提出了矩阵理论第二大猜想：给定一个封闭的共轭n元复数组 $\left\{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\right\}$，能否找出使得此n元复数组是某个非负矩阵谱的充分必要条件。文献 [4] [7] [8] [9] [10] 研究了建立当矩阵的阶数小于等于5阶时，有解的充分必要条件；当阶数高于5阶时，文献 [9] [10] 研究了有解的充分条件。其中Perron-Frobenius定理的运用，使得这个问题的解决取得了很大的进步；关于解的存在性比较完整的研究结果可以参考M.T. Chu和G.H. Golub合著的文献 [11]。

本文关注的重点是NIEP的数值算法。关于已知部分特征对的NIEP，白正简老师等人在 [2] 中提出了非光滑牛顿法，并证明了算法局部收敛，且算法收敛时收敛速度是二阶的。进一步，白正简老师等人在文献 [1] 中将这类问题转化为黎曼流形上的优化问题，提出了共轭梯度算法，并证明了算法的收敛性。最近，杨丹在 [12] 中，通过将这类问题转化为凸可行性问题，提出了交替投影算法，并证明了算法线性收敛。关于经典的NIEP，即已知矩阵的全体特征值，求非负矩阵，Oris等人在文献 [8] 中将非负逆特征值问题转化为求解两个集合交集的可行性问题，从而提出了交替投影算法，并在一定条件下建立了算法的收敛性。受这一思想的启发，党婵娟在文献 [13] 中研究了StNIEP。通过将该问题转化为求解可行性问题，提出了交替投影算法和Average交替投影算法，并借助黎曼流形上投影算子的性质研究了算法的收敛性和收敛速度。值得一提的是，Chu等人在文献 [7] 将SNIEP转化为黎曼流形上的优化问题，通过建立相应变量的微分方程，提出了“梯度流”算法。该问题本质上是将SNIEP转化为微分方程来求解，当文中的 $\omega$ 极限点满足一定条件时，算法收敛(详见 [7])。最近，Chu等人在 [11] 运用类似的技术研究了一般NIEP的梯度流算法，同时运用到求解StNIEP。

自从1972年Luenberger提出黎曼流形上沿测地线的梯度下降算法 [14]，黎曼流形上的优化算法的研究受到了很大的关注。其他黎曼流形上优化算法的研究，包括更一般的梯度算法，信赖域算法，牛顿法，次梯度法等见 [15] [16] - [21] 及其参考文献。正如这些文献中指出，运用黎曼流形技术求解欧氏空间的优化问题，不论在理论上还是实际应用中都有重要的意义。其主要优点包括可将欧氏空间中的约束优化问题转化为无约束优化问题；将病态问题转化为非病态问题；将非凸优化问题转化为凸优化问题。本文在第三节考虑的优化问题就运用把欧氏空间中的约束优化问题转化为无约束优化问题这一优点。受到Chu等人在文献 [7] 研究的启发，基于该文献提出的黎曼优化模型，我们分别提出黎曼梯度算法和黎曼共轭梯度算法求解SNIEP。通过对SNIEP的分析，分别建立了两类算法的全局收敛性(算法的收敛性不依赖于初始点的选取)。正如我们上面指出，已有求解SNIEP的算法，均不具备全局收敛性。文献 [7] 中的梯度流算法当 $\omega$ 极限点满足一定条件时收敛，而文献 [8] 中交替投影算法在初始点满足一定条件下收敛。最后通过数值例子，我们验证了算法的收敛性。在数值实验中，将Stiefel流形上的拉回映射和平移映射应用到算法中，解决了乘积黎曼流形上的拉回映射和平移映射的计算问题。通过数值例子可以看出，求解高阶的SNIEP时，黎曼共轭梯度算法的收敛效率高于黎曼梯度算法。

本文结构如下。在下一节，我们将介绍一些黎曼流形上的相关知识。第三节分别提出来求解SNIEP的黎曼梯度下降算法和黎曼共轭梯度算法。第四节，分析得到了算法的全局收敛性。最后，通过第五节的数值实验，验证并比较了两种算法收敛性。

2. 预备知识

这节主要介绍黎曼流形上的一些基本的概念和结论，它们可以在黎曼流形的书上找到，例如 [15] [22] [23]。

设 $\mathcal{M}$ 是一个Hausdorff拓扑空间，若 $\mathcal{M}$ 的每一个点x都有一个开邻域 $U\subset \mathcal{M}$，使得U和n维欧氏空间 ${ℝ}^n$ 中的一个开子集是同胚的，则称 $\mathcal{M}$ 是一个n维拓扑流形。若 $\mathcal{M}$ 上有 $C^{\infty }$ 微分结构，则称 $\mathcal{M}$ 为 $C^{\infty }$ 微分流形(或光滑流形)。以下我们总假设 $\mathcal{M}$ 是光滑流形。设 $x\in \mathcal{M}$，则在点x处的所有切向量组成的集合称为流形 $\mathcal{M}$ 在点x处的切空间，记为 $T_x\mathcal{M}$，而 $T\mathcal{M}:=\underset{x\in \mathcal{M}}{\cup }{T}_x\mathcal{M}$ 称为切丛。设g是 $\mathcal{M}$ 上的一个光滑的二阶协变张量场。如果g是对称、正定的，则称g为 $\mathcal{M}$ 上的黎曼度量。若在 $\mathcal{M}$ 上指定一个黎曼度量g，则称 $\left(\mathcal{M},g\right)$ 为一个黎曼流形，简记为 $\mathcal{M}$。本文也用符号 $\langle \cdot ,\cdot \rangle$ 来表示在光滑流形上指定的黎曼度量g。在 $\mathcal{M}$ 上存在唯一一个与度量g相容的黎曼联络D。设 $I\subset ℝ$ 是 $ℝ$ 中一个区间， $\gamma :I\to \mathcal{M}$ 为 $\mathcal{M}$ 上一条光滑

曲线，如果满足 $\frac{D}{\text{d}t}\left(\frac{\text{d}\gamma }{\text{d}t}\right)=0$，则称其为 $\mathcal{M}$ 上一条测地线。设 $\mathcal{M}$ 是完备的黎曼流形，则对任意 $x\in \mathcal{M}$，任意 $\xi \in T_x\mathcal{M}$，则存在唯一测地线 $\gamma =\gamma \left(t;x,\xi \right):ℝ\to \mathcal{M}$，满足：

${\gamma \left(0\right)=x,\frac{\text{d}\gamma }{\text{d}t}|}_{t=0}=\xi$。

设 $x\in \mathcal{M}$，定义点x处的指数映射为：

。

我们知道欧氏空间 ${ℝ}^n$ 作为特殊的黎曼流形，其每一点处的切空间都是它自己。 ${ℝ}^n$ 中点和点以及点和切向量之间可以直接进行线性运算。但是一般的黎曼流形上不存在这样的线性运算，为进行类似的运算，我们需要引入拉回映射、平移映射作等。以下拉回映射的概念是指数映射的推广，它将切空间 $T_x\mathcal{M}$ 上的点一一映射到流形 $\mathcal{M}$ 上，这一概念最先由Adler [24] 提出。

定义2.1称流形 $\mathcal{M}$ 上的光滑映射是拉回映射，如果 $\mathcal{R}$ 在 $T_x\mathcal{M}$ 上的映射 ${\mathcal{R}}_x$ 满足：

1) ${\mathcal{R}}_x\left(0_x\right)=x$， $0_x$ 是 $T_x\mathcal{M}$ 的零元；

2) 若，则 ${\mathcal{R}}_x$ 满足 $D{\mathcal{R}}_x\left(0_x\right)={\text{id}}_{T_x\mathcal{M}}$， ${\text{id}}_{T_x\mathcal{M}}$ 是 $T_x\mathcal{M}$ 上的恒等映射。

可以验证，指数映射是拉回映射。除指数映射以外，还有其他拉回映射，例如在球面上，， $\mathcal{R}$ 在 $T_x\mathcal{M}$ 上的映射 ${\mathcal{R}}_x$ 定义为：

，

则 $\mathcal{R}$ 是球面上的拉回映射。

为在黎曼流形 $\mathcal{M}$ 上不同点的切空间之间实现加(减)法运算，定义平移映射 [15]。为此，记 $T\mathcal{M}\oplus T\mathcal{M}=\left\{\left(\eta ,\xi \right):\eta ,\xi \in T_x\mathcal{M},x\in \mathcal{M}\right\}$。

定义2.2称黎曼流形 $\mathcal{M}$ 上的光滑映射是 $\mathcal{M}$ 上的平移映射，如果 $\mathcal{T}$ 满足：

1) 对于任意 $\xi ,\eta \in T_x\mathcal{M},x\in \mathcal{M}$，存在一个拉回映射 $\mathcal{R}$，使得 ${\mathcal{T}}_{\eta }\left(\xi \right)\in T_{{\mathcal{R}}_x\left(\eta \right)}\mathcal{M}$ ；

2) 对于任意 $\xi \in T_x\mathcal{M},x\in \mathcal{M}$，都有 ${\mathcal{T}}_{0_x}\left(\xi \right)=\xi$ ；

3) 对于任意 $\xi ,\zeta \in T_x\mathcal{M},x\in \mathcal{M}$，都有 ${\mathcal{T}}_{\eta }\left(a\xi +b\zeta \right)=a{\mathcal{T}}_{\eta }\left(\xi \right)+b{\mathcal{T}}_{\eta }\left(\zeta \right)$。

特别地，我们介绍Stiefel流形 $\mathcal{M}=S\text{t}\left(p,n\right)=\left\{X\in {ℝ}^{n\times p}:X^{\text{T}}X=I_p\right\}\left(p\ge 2\right)$ 及其上的拉回映射和平移映射等，这些结论均可以在 [15] 中找到。设，在点X处切空间为

$T_X\mathcal{M}=\left\{Z\in {ℝ}^{n\times p}:X^{\text{T}}Z+Z^{\text{T}}X=0\right\}=\left\{X\Omega +X_{\perp }K:{\Omega }^{\text{T}}=-\Omega ,K\in {ℝ}^{\left(n-p\right)\times p}\right\}$，

黎曼度量为

，

其中 $\mathrm{trace}\left(A)\right)$ 表示矩阵A的迹。设。在点X处指数映射为：

${\text{Exp}}_X\left(t\right)=\left[X,Z\right]\exp \left(t\left[\begin{array}{cc}A\left(0\right)& -S\left(0\right)\\ I_p& A\left(0\right)\end{array}\right]\right)\left[\begin{array}{c}{I}_p\\ 0\end{array}\right]\exp \left(-A\left(0\right)t\right)$，

其中 $A\left(0\right):=X^{\text{T}}Z,S\left(0\right)=Z^{\text{T}}Z$。基于极分解的拉回映射定义为：

${\mathcal{R}}_X\left(Z\right)=\left(X+Z\right){\left(I_p+Z^{\text{T}}Z\right)}^{-\frac{1}{2}}$ ；

基于QR分解的拉回映射可定义为：

${\mathcal{R}}_X\left(Z\right)=\text{qf}\left(X+Z\right)$， (2.1)

其中 $\text{qf}\left(X+Z\right)$ 是对矩阵 $X+Z$ 进行QR分解后所得到的正交矩阵Q [15]。在流形上可以沿着测地线进行线搜索，也可以直接使用相应的拉回映射，但它们在不同的流形上的计算效率是不一样的。在球面上无论是沿测地线进行线搜索还是使用拉回映射，计算效率相差不大 [3] [25]。一般而言，即使同一个流形上，不同的拉回映射也可能会有不同的计算效率。设， $Z_1,Z_2\in T_X\mathcal{M}$。在后面的数值实验中，我们采用以下Stiefel流形上的平移映射 [22]：

${\mathcal{T}}_{Z_2}{Z}_1=\left(I-Y{Y}^{\text{T}}\right)Z_1+Y\text{skew}\left(Y^{\text{T}}{Z}_1\right)\in T_Y\mathcal{M}$， (2.2)

其中 $Y={\mathcal{R}}_X\left(Z_2\right),\text{skew}\left(Y^{\text{T}}{Z}_1\right):=\frac{1}{2}\left(Y^{\text{T}}{Z}_1-Z_1^{\text{T}}Y\right)$。注意，其中拉回映射 $\mathcal{R}$ 可以是上面指出的拉回映射中的任意一个。

设 $f:\mathcal{M}\to ℝ$ 是黎曼流形 $\mathcal{M}$ 上的光滑函数，定义f的梯度为一个光滑切向量场 $\text{grad}f$ 如下：

$g\left(\text{grad}f\left(x\right),\xi \right)=\langle \text{grad}f\left(x\right),\xi \rangle =\xi \left(f\right)\forall x\in \mathcal{M},\xi \in T_x\mathcal{M}$。

设黎曼流形 $\mathcal{M}$ 是黎曼流形 $\bar{\mathcal{M}}$ 的子流形， $x\in \mathcal{M}$， $\bar{f}:\bar{\mathcal{M}}\to ℝ$ 是 $\bar{\mathcal{M}}$ 的光滑函数。设 $f:\mathcal{M}\to ℝ$ 是 $\bar{f}$ 在 $\mathcal{M}$ 上的限制，则f在 $\mathcal{M}$ 上梯度函数等于 $\bar{f}$ 在 $\bar{\mathcal{M}}$ 上的梯度在 $T_x\mathcal{M}$ 上的投影，即有：

$\text{grad}f\left(x\right)=P_x\text{grad}\bar{f}\left(x\right)$，

其中 $P_x$ 是从切空间 $T_x\bar{\mathcal{M}}$ 到切空间的投影算子。本文中的范数“ $\Vert \cdot \Vert$ ”均表示为Frobenius矩阵范数。

3. 黎曼梯度下降算法及黎曼共轭梯度算法

考虑黎曼流形上的优化问题：

$\underset{X\in \mathcal{M}}{\min }f\left(X\right)$，

其中 $\mathcal{M}$ 是一个黎曼流形， $f:\mathcal{M}\to ℝ$ 是定义在 $\mathcal{M}$ 上的光滑函数。正如我们在引言部分指出，黎曼流形上的一阶优化算法的研究取得了丰富的研究结果。特别地，文献 [1] [15] 分别研究了黎曼流形上带Armijo步长的梯度下降算法和共轭梯度算法，具体见下面的算法1和算法2。

算法1. 采用Armijo步长的黎曼梯度下降算法

算法2. 采用Armijo步长的黎曼共轭梯度算法

关于这两个算法的收敛性见下面的命题3.1和3.2，它们分别来自文献 [15] [18] [19] [25] 和 [1] [19]。

命题3.1假设 $f:\mathcal{M}\to ℝ$ 是一个光滑映射， $\left\{X_k\right\}$ 是在f作用下由算法1产生的无穷点列，且满足：

① 对于任意给定点 $X_0\in \mathcal{M}$， $\Omega =\left\{X\in \mathcal{M}:f\left(X\right)\le f\left(X_0\right)\right\}$ 是紧集；

② 梯度函数 $\text{grad}f:T_X\mathcal{M}\to T_X\mathcal{M}$ 在 $0_X\in T_X\mathcal{M}$ 处Lipschitz连续，即存在一个常数

$L>0,X_1,X_2\in \mathcal{M}$，使得：

$\Vert \text{grad}f\left(X_1\right)-\text{grad}f\left(X_2\right)\Vert \le L\cdot \text{dist}\left(X_1,X_2\right)$，

其中 $\text{dist}\left(X_1,X_2\right)$ 是 $X_1,X_2$ 两点的黎曼距离。则有

$\underset{k\to \infty }{\lim \inf }\Vert \text{grad}f\left(X_k\right)\Vert =0$ (3.1)

命题3.2假设 $f:\mathcal{M}\to ℝ$ 是一个光滑映射， $\left\{X_k\right\}$ 是在f作用下由算法2产生的无穷点列，满足命题3.1中的条件①和②，且算法2中的参数还满足条件

③ $1-\frac{{\mu }_1}{{\mu }_2}>\frac{2{\delta }_1}{1+{\delta }_1}$，且 $\frac{{\mu }_1}{{\mu }_3}<1$， (3.2)

则有(3.1)成立。

4. 对称非负逆特征值问题

本文研究求解SNIEP的算法，即求一个对称非负矩阵，使得它的特征值是给定的n个常数 $\left\{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\right\}$。我们将通过SNIEP转化为将黎曼流形上优化问题，并提出黎曼梯度算法和黎曼共轭梯度算法进行求解。为此，先引进一些记号。记 $\Lambda :=\text{diag}\left\{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\right\}$。记所有的n阶正交矩阵构成的Stiefel流形为 $\mathcal{O}\left(n\right)$ 集合，即

$\mathcal{O}\left(n\right):=\left\{Q\in {ℝ}^{n\times n}:Q^{\text{T}}Q=I_n\right\}$。

设 $X\in \mathcal{O}\left(n\right)$，则流形 $\mathcal{O}\left(n\right)$ 在点X处的切空间是

$T_X\mathcal{O}\left(n\right)=\left\{Z\in {ℝ}^{n\times n}:X^{\text{T}}Z+Z^{\text{T}}X=0\right\}$。

黎曼度量为

$g\left(Z_1,Z_2\right)=\langle Z_1,Z_2\rangle =\mathrm{trace}\left(Z_1^{\text{T}}{Z}_2\right)\forall Z_1,Z_2\in T_X\mathcal{O}\left(n\right)$。

记n阶实对称矩阵组成的集合 $\mathcal{S}\left(n\right):=\left\{B\in {ℝ}^{n\times n}:B=B^{\text{T}}\right\}$，显然 $\mathcal{S}\left(n\right)$ 是欧氏空间，也是特殊的黎曼流形。设 $X\in \mathcal{S}\left(n\right)$，则 $\mathcal{S}\left(n\right)$ 在X点处的切空间为

$T_X\mathcal{S}\left(n\right)=\mathcal{S}\left(n\right)$。

在点X处的黎曼度量为

$g\left(W_1,W_2\right)=\langle W_1,W_2\rangle =\mathrm{trace}\left(W_1^{\text{T}}{W}_2\right)\forall W_1,W_2\in T_X\mathcal{S}\left(n\right)$。

显然乘积空间 $\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 也构成一个黎曼流形，在任意一点 $\left(Q,R\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 处的切空间为

$T_{\left(Q,R\right)}\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)=T_Q\mathcal{O}\left(n\right)\times T_R\mathcal{S}\left(n\right)=\left(T_Q\mathcal{O}\left(n\right),T_R\mathcal{S}\left(n\right)\right)$。

设 $\left(Q,R\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$。在点 $\left(Q,R\right)$ 的黎曼度量为：对任意 $\left(Z_1,W_1\right),\left(Z_2,W_2\right)\in T_{\left(Q,R\right)}\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$，

$g_{\left(Q,R\right)}\left(\left(Z_1,W_1\right),\left(Z_2,W_2\right)\right)=\langle \left(Z_1,W_1\right),\left(Z_2,W_2\right)\rangle =\text{trace}\left(Z_1^{\text{T}}{Z}_2\right)+\text{trace}\left(W_1^{\text{T}}{W}_2\right)$。

对于黎曼流形 $\mathcal{O}\left(n\right)$，本文采用了基于QR分解的拉回映射(2.1)和平移映射(2.2)。设 $\left(Q,R\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$， $\left(Z_1,W_1\right),\left(Z_2,W_2\right)\in T_{\left(Q,R\right)}\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$，则可以分别定义流形 $\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 上的拉回映射和平移映射如下：

${\mathcal{R}}_{\left(Q,R\right)}\left(Z_1,W_1\right)=\left(\text{qf}\left(Q+Z_1\right),R+W_1\right)$，

${\mathcal{T}}_{\left(Z_1,W_1\right)}\left(Z_2,W_2\right)=\left(\frac{1}{2}\left(Z_1-Y{Z}_1^{\text{T}}Y\right),W_2\right)$，

其中 $Y:={\mathcal{R}}_Q\left(Z_2\right)$。

文献 [1] 将SNIEP转化成了下面的黎曼优化问题：

$\begin{array}{l}\min f\left(Q,R\right):=\frac{1}{2}{\Vert Q^{\text{T}}\Lambda Q-R\ast R\Vert }^2\\ \text{s}\text{.t}\text{.}\left(Q,R\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)\end{array}$ (4.1)

其中 $\Vert \cdot \Vert$ 表示矩阵的F范数， $R\ast R$ 表示矩阵R和R的Hadamard乘积(即 $R\ast R=\left(R_{ij}^2\right)$， $R=\left(R_{ij}\right)$ )。当SNIEP有解时，SNIEP等价于黎曼流形 $\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 上优化问题(4.1)。同时文献 [7] 给出了目标函数f在流形 $\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 上的梯度表达式如下：对于任意的 $\left(Q,R\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$，都有

$\text{grad}f\left(Q,R\right)=\left(Q\left(R\ast R\right)Q^{\text{T}}\Lambda Q-\Lambda Q\left(R\ast R\right),-2\left(Q^{\text{T}}\Lambda Q-R\ast R\right)\ast R\right)$。 (4.2)

Chu等人在 [7] 提出了用“梯度流”的方法求解问题(4.1)。具体地，将 $X\left(t\right)=Q{\left(t\right)}^{\text{T}}\Lambda Q\left(t\right)$ 、 $Y\left(t\right)=R\left(t\right)\ast R\left(t\right)$ 看作关于t的“流”，通过建立 $Q\left(t\right)$ 和 $R\left(t\right)$ 满足的微分方程，运用微分方程的求解方法求解其极限点。如果轨道 $\left(X\left(t\right),Y\left(t\right)\right)$ 的 $\omega$ -极限集包含一个类型为 $\left(\hat{X},\hat{X}\right)$ 的单点，则 $\left(X\left(t\right),Y\left(t\right)\right)$ 收敛到 $\left(\hat{X},\hat{X}\right)$。容易看出，“梯度流”的求解思路本质上是最速下降法的思路，但不是所有梯度“流”都会收敛到问题(4.1)的解。基于黎曼流形上的优化模型(4.1)，我们提出采用线搜索的黎曼优化算法1和算法2进行求解。以下是这两个算法运用到求解问题(4.1)的收敛性分析，为此我们先证明下面的引理。

引理4.1在模型(4.1)中，对于任意给定点 $\left(Q_0,R_0\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$，以下水平集是紧集：

$\Omega =\left\{\left(Q,R\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right):f\left(Q,R\right)\le f\left(Q_0,R_0\right)\right\}$。

证明：设 $\left(Q,R\right)\in \Omega$，则

$\Vert Q\Vert =\sqrt{\text{trace}\left(Q^{\text{T}}Q\right)}=\sqrt{n}$， (4.3)

$\Vert Q^{\text{T}}\Lambda Q-R\ast R\Vert \le \sqrt{2f\left(Q_0,R_0\right)}$。 (4.4)

由柯西不等式，

$\Vert R\Vert =\sqrt{{\displaystyle \underset{i,j}{\sum }{r}_{ij}^2}}=\sqrt{{\displaystyle \underset{i,j}{\sum }\left(1\cdot r_{ij}^2\right)}}\le \sqrt{n^2\cdot {\displaystyle \underset{i,j}{\sum }{r}_{ij}^4}}=n\cdot \Vert R\ast R\Vert$。 (4.5)

又由矩阵F范数的酉不变性 [26] 可得

$\Vert Q^{\text{T}}\Lambda Q\Vert =\Vert \Lambda \Vert$。

所以由(4.4)，得

$\Vert R\ast R\Vert \le \Vert Q^{\text{T}}\Lambda Q-R\ast R\Vert +\Vert Q^{\text{T}}\Lambda Q\Vert \le \sqrt{2f\left(Q_0,R_0\right)}+\Vert \Lambda \Vert$，

所以由(4.5)

$\Vert R\Vert \le n\cdot \left(\sqrt{2f\left(Q_0,R_0\right)}+\Vert \Lambda \Vert \right)$，

于是由(4.3)及上式可得

$\Vert \left(Q,R\right)\Vert =\Vert Q\Vert +\Vert R\Vert \le n\cdot \left(\sqrt{2f\left(Q_0,R_0\right)}+\Vert \Lambda \Vert \right)+\sqrt{n}$，

故水平集 $\Omega$ 是有界的。又f在 $\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 上显然是连续的，故水平集 $\Omega$ 是闭集，从而 $\Omega$ 是紧集。

引理4.2设 $\left(Q_0,R_0\right)\in \mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$，则 $\text{grad}f$ 在水平集 $\Omega$ 上是Lipschitz连续的( $\Omega$ 如引理4.1中定义)。

证明：由目标函数f的定义，可知f是定义在 $\mathcal{O}\left(n\right)\times \mathcal{S}\left(n\right)$ 上的光滑函数。又由引理4.1知 $\Omega$ 是紧集，f的Hessian矩阵在 $\Omega$ 上有界。再由微分中值定理， $\text{grad}f$ 在 $\Omega$ 上是Lipschitz连续。

由命题2.1，2.2和引理4.1，4.2，可得以下定理。

定理4.1设 $\left\{X_k\right\}$ 是运用算法1或者算法2求解问题(4.1)产生的无穷点列，且算法2中参数满足命题2.2中条件③，则有

$\underset{k\to \infty }{\lim \inf }\Vert \text{grad}f\left(X_k\right)\Vert =0$。

5. 数值实验

本节利用数值实验来验证算法1和算法2有效性，并比较这两种算法之间的收敛性。停机准则： ${\Vert \text{grad}f\left(Q,R\right)\Vert }^2\le {10}^{-12}$，实验所用软件是Matlab2018a。实验中初始矩阵 $Q_0=I_n$ 是n阶单位矩阵， $R_0=ones\left(n\right)$ 是所有元素都是1的n阶矩阵。实验中的参数取值情况如下： $\alpha =\rho =0.5$， $\delta ={\delta }_1=0.1$， ${\delta }_2=10$。为了满足命题3.2的条件③ (即(3.2)式)，取 ${\mu }_1=1,{\mu }_2=2,{\mu }_3=3$。下图1~4所示的是不同阶数，两种算法的收敛比较(横坐标是迭代次数，纵坐标是函数值)。

Figure 1. n = 4

图1. n = 4

Figure 2. n = 6

图2. n = 6

Figure 3. n = 20

图3. n = 20

Figure 4. n = 100

图4. n = 100

为了进一步比较黎曼梯度下降算法和黎曼共轭梯度算法的收敛效率，在选取和以上实验相同参数的情况下，将函数值的精确度Res取为1.0E−13；矩阵的阶数n依次取为6，10，20，50，100，迭代时间都是运行程序10次后取的平均值。其中i是迭代的次数，T是迭代时间(单位是秒)，Res表示函数值f的精确度，n是矩阵的阶数。见表1。

Table 1. The experimental results of Algorithm 1 and Algorithm 2 at different orders

表1. 算法1和算法2在不同阶数下的实验结果

从以上各图可以看出，在相同的迭代次数时，黎曼梯度下降算法比的函数值下降量更多。但是从上表可以看出，在阶数较高的情形下( $n\ge 50$ )，黎曼共轭梯度算法达到相同精度的运行时间比黎曼梯度下降算法更短。这表明，在求解高阶SNIEP时，黎曼共轭梯度算法比黎曼梯度下降算法更有效。

6. 结论

本文提出了用黎曼梯度下降算法(算法1)和黎曼共轭梯度算法(算法2)求解对称非负逆特征值问题，并建立了算法的收敛性。在数值实验中比较了这两种算法的收敛效率，实验的结果是：在高阶的情况下，黎曼共轭梯度算法比黎曼梯度下降算法收敛效率更高。接下来我们将研究用黎曼二阶算法求解对称非负逆特征值问题，并比较各种黎曼算法的收敛效率。

基金项目

本文受国家自然科学基金(12161017)和贵州省省级科技计划项目(ZK[2022]110)资助项目资助。

NOTES

^*通讯作者。