交替投影算法求解非负逆特征值问题

Alternating Projection Method for Solving Nonnegative Inverse Eigenvalue Problems

摘要: 通过把给定部分特征对的非负逆特征值问题转化为一个凸可行性问题，提出交替投影算法求解该问题。建立了这一算法的线性收敛性。最后，通过数值例子，比较了交替投影算法和非光滑牛顿法(白等人2011年提出)的收敛效率。数值实验结果表明，交替投影算法总是能收敛到问题的解，而非光滑牛顿法在一些情形下求不出解。此外，交替投影算法收敛的效率也比非光滑牛顿法高。

Abstract: The alternating projection method is proposed for solving nonnegative inverse eigenvalue prob-lems with partial eigendata, by reformulating the problem as a convex feasibility problem. The (linear) convergnece property of the method is established. At last, some numerical experiments are provided to compare the method with the non-smooth Newton algorithm proposed by Bai et al. in 2011. Numerical experiment results show that the alternative projection algorithm always converges to the solution of the problem, while the non-smooth Newton method fails to find the solution in some cases. In addition, the alternative projection algorithm is more efficient than the non-smooth Newton method.

1. 引言

非负逆特征值问题来源非常广泛，它主要来自于离散的数学物理反问题、系统参数识别、地震断层成像技术、主成分分析与勘测、结构分析、电路理论、机械系统模拟等许多应用领域，有着重要的应用背景 [1] - [6] 。这类问题提出来的几十年里，吸引了大量学者，取得到了重要研究成果。非负逆特征值问题可以分为以下两大类：1) 已知n个特征值 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$ ，求一个非负n阶非负矩阵¹ (或者非负对称矩阵/

随机矩阵)；2) 已知部分特征对 ${\left\{{\lambda }_k,x_k\right\}}_{k=1}^p\left(p\le n\right)$ ，其中 ${\lambda }_k$ 是特征值， $x_k$ 是其对应的特征向量，求一个n

阶非负矩阵(或者非负对称矩阵/随机矩阵)。关于这两类问题的研究主要包括两方面 [7] ，一方面是关于问题解的存在性；另一反面是求解这些问题的数值算法。例如，求解第一类问题的等光谱流动算法 [8] 以及交替投影算法 [9] 等。关于第二类问题的数值算法求解见文献 [10] [11] [12] ，其中文献 [10] 通过把非负逆特征值问题转化为优化问题，提出了非光滑牛顿法。在满足某种非退化的条件下，证明了算法是二阶收敛的。但该算法在求解子问题的计算量大，收敛的效率并不理想。

本文主要研究第二类非负逆特征值问题。通过把这一问题转化为凸可行性问题，提出了交替投影算法，并给出了投影的计算式。同时，我们证明了算法的(线性)收敛性。最后，通过数值例子比较了交替投影算法和非光滑牛顿法的收敛效率。数值实验结果表明，交替投影算法总能收敛到问题的解，而非光滑牛顿法在一些情形下不收敛。此外，非光滑牛顿法虽然在一定条件下二阶收敛，但是因为子问题计算量大，非光滑牛顿法的实际运算效率不如交替投影算法。

本文的结构如下。第二节是预备知识，包括一些相关的概念和记号。文章的主要结果在第三节，我们提出求解非负逆特征问题的交替投影算法，并证明了算法的线性收敛性。第四节是数值实验及与非光滑牛顿法的比较。

2. 预备知识

这节主要介绍矩阵分析中的一些符号，概念和结论，参考 [13] [14] [15] 。设 $R^n$ 和 $C^n$ 分别表示n维实向量空间和复向量空间， $R^{m\times n}$ 和 $C^{m\times n}$ 分别表示 $m\times n$ 阶实矩阵空间和 $m\times n$ 阶复矩阵空间。矩阵 $A\in R^{m\times n}$ 的

转置，共轭转置，迹分别记为 $A^{\text{T}}$ ， $A^H$ ， $\mathrm{tr}\left(A\right)$ ，其中 $\mathrm{tr}\left(A\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{a}_{ii}}$ 。在 $R^{n\times n}$ 中，定义内积 $\langle \cdot ,\cdot \rangle$ 如下：

$\langle A,B\rangle =\mathrm{tr}\left(A^{\text{T}}B\right),\forall A,B\in R^{n\times n}.$

其诱导Frobenius模记为 ${\Vert \cdot \Vert }_F$ ，即

${\Vert A\Vert }_F=\sqrt{\langle A,A\rangle }=\sqrt{{\displaystyle \underset{i,j=1}{\overset{n}{\sum }}{a}_{ij}^2}},A\in R^{n\times n}.$

设矩阵 $A\in C^{m\times n}$ 。如果矩阵 $X\in C^{m\times n}$ 满足以下四个方程：

$AXA=A,XAX=X,{\left(AX\right)}^H=AX,{\left(XA\right)}^H=XA.$

则称矩阵X为矩阵A的广义逆矩阵，记作 $A^{+}$ 。

假设 $S\subseteq R^{n\times n}$ 为 $R^{n\times n}$ 中非空子集。设 $X\in R^{n\times n}$ ，点X到集合S上的距离和投影分别记为 $d_S(\cdot ),P_S(\cdot )$ 即

$d_S\left(X\right):=\underset{Y\in S}{\inf }\Vert X-Y\Vert$

和

$P_S\left(X\right)=\underset{Y\in S}{\arg \min }\Vert X-Y\Vert =\left\{Y\in S|d_S\left(X\right):=\Vert X-Y\Vert \right\}.$

如果任意两点 $x_1,x_2\in S$ ，都有

${\lambda }_1{x}_1+\left(1-\lambda \right)x_2\in S,\forall \lambda \in \left(0,1\right),$

则称S为 $R^{n\times n}$ 中的凸集。设 $A\in R^{n\times n},b\in R$ ，我们称集合 $H=\left\{X\in R^{n\times n}|\langle A,X\rangle \le b\right\}$ 为 $R^{n\times n}$ 中的半空间，几个半空间的交称为 $R^{n\times n}$ 中的多面体，容易验证半空间和多面体都是凸集。

3. 非负逆特征值问题及交替投影算法

本文主要研究以下非负逆特征值问题(下面简记为NIEP)。已知p个特征对 ${\left\{{\lambda }_k,x_k\right\}}_{k=1}^p\left(p\le n\right)$ 是A的p

个特征对，其中 ${\lambda }_k\in C$ ， $x_k\in C^n$ 是特征值 ${\lambda }_k$ 对应的特征向量，求一个非负矩阵A。为此我们先引进一些记号。不失一般性，我们不妨假设前2s个是共轭特征对，其余是实特征对，即对任意 $1\le i\le s$ ，

${\lambda }_{2i-1}=a_i+b_i\sqrt{-1},{\lambda }_{2i}=a_i-b_i\sqrt{-1},x_{2i-1}=x_{iR}+x_{iI}\sqrt{-1},x_{2i}=x_{iR}-x_{iI}\sqrt{-1},$

其中 $a_i,b_i\in R,x_{iR},x_{iI}\in R^n$ 。记

$\Lambda =diag\left({\lambda }_1^{\left[2\right]},\cdots ,{\lambda }_s^{\left[2\right]},{\lambda }_{2s+1},\cdots ,{\lambda }_p\right)\in R^{P\times P}$

和

$X=\left[x_{1R},x_{1I},\cdots ,x_{sR},x_{sI},x_{2s+1},\cdots ,x_p\right]\in R^{n\times p},$

其中 ${\lambda }_i^{\left[2\right]}=\left[\begin{array}{cc}{a}_i& b_i\\ -b_i& a_i\end{array}\right]$ 。则NIEP是求一个非负矩阵 $A\in R^{n\times n}$ ，使得

$AX=X\Lambda .$ (1)

设

$C_1=\left\{A\in R^{n\times n}|AX=X\Lambda \right\},C_2=R_{+}^{n\times n}$ (2)

( $R_{+}^{n\times n}$ 是全体n阶非负矩阵的集合)。容易验证 $C_1,C_2$ 都是多面体，从而是凸集.则NIEP可以转化为以下凸可行性问题：

求矩阵 $A\in R^{n\times n}$ ，使得 $A\in C_1\cap C_2$ 。 (3)

通过把NIEP转化为问题(3)，我们用以下交替投影算法求解NIEP。关于交替投影算法的其他应用和研究结果可见参考文献 [16] 。

算法3.1 (交替投影算法)

Step 0 给定初始矩阵 $A^0\in R^{n\times n}$ ， $\epsilon >0$ ，令 $k=0$ 。

Step 1令 $A^{k+1}=P_{C_2}\left(P_{C_1}\left(Y^k\right)\right)$ 。

Step 2 若 ${\Vert A^{k+1}X-X\Lambda \Vert }_F\le \epsilon$ ，算法终止；否则，进入下一步。

Step 3 令 $k=k+1$ ，转Step 1。

下面分别给出 $P_{C_1}(\cdot )$ 和 $P_{C_2}(\cdot )$ 的计算式。由定义容易验证，对任意 $A=\left(a_{ij}\right)\in R^{n\times n}$ ，有 $P_{C_2}\left(A\right)={\left({\bar{a}}_{ij}\right)}_{n\times n}$ ，其中

${\bar{a}}_{ij}=\{\begin{array}{l}{a}_{ij},a_{ij}\ge 0\\ 0,a_{ij}<0\end{array},i,j=1,2,\cdots ,n.$

为得到 $P_{C_1}(\cdot )$ 的计算式，我们引用以下结论，见 [10] 。

引理3.1设 $B_1\in R^{q\times m},B_2\in R^{n\times t},B_3\in R^{q\times t}$ ，其中 $p,q,t$ 为正整数。则集合 $C=\left\{A\in R^{m\times n}|B_{1}A{B}_2=B_3\right\}$ 非空的充分必要条件是 $B_1,B_2,B_3$ 满足 $B_1{B}_1^{+}{B}_3{B}_2^{+}{B}_2=B_3$ 。并且当C非空时有

$P_C\left(Y\right)=B_1^{+}{B}_3{B}_2^{+}+A-B_1^{+}{B}_{1}A{B}_2{B}_2^{+},\forall Y\in R^{m\times n}.$

引理3.2设集合 $C_1$ 由(2)定义，则 $C_1$ 非空的充分必要条件是 $X\Lambda X^{+}X=X\Lambda$ 。并且当 $C_1$ 非空时有

$P_{C_1}\left(Y\right)=X\Lambda X^{+}+Y\left(I-X{X}^{+}\right),\forall Y\in R^{n\times n}.$

证：由定义有 $C_1=\left\{A\in R^{n\times n}|EAX=X\Lambda \right\}$ ，其中E是n阶单位阵。令 $B_1=E,B_2=X,B_3=X\Lambda$ ，由引理3.1， $C_1$ 非空的充分必要条件是 $B_1{B}_1^{+}{B}_3{B}_2^{+}{B}_2=B_3$ ，即 $X\Lambda X^{+}X=X\Lambda$ 。并且此时对任意 $Y\in R^{n\times n}$ ，有 $P_{C_1}\left(Y\right)=B_1^{+}{B}_3{B}_2^{+}+Y-B_1^{+}{B}_{1}Y{B}_2{B}_2^{+}=X\Lambda X^{+}+Y\left(I-X{X}^{+}\right)$ 。证毕。

定理3.3设NIEP解集非空， $\left\{A^k\right\}$ 是算法3.1产生的点列。则 $\left\{A^k\right\}$ 线性收敛到NIEP的一个的解。

证：用定义容易验证 $C_1,C_2$ 是多面体。从而由 [17] (或 [16] )，知 $\left\{C_1,C_2\right\}$ 线性正则。再由 [17] (或 [16] )，可知算法3.1产生的点列 $\left\{A^k\right\}$ 线性收敛到问题(3)的一个解。从而 $\left\{A^k\right\}$ 线性收敛到NIEP的一个解。证毕。

4. 数值实验结果

在这一节，我们将用数值实验验证算法3.1的收敛性，并比较该算法和非光滑牛顿法 [10] 的收敛效率。注意到，在文献 [10] 中，作者证明了非光滑牛顿法在某种非退化的条件下二阶收敛。所以非光滑牛顿法不能保证对所有初始点算法收敛。此外，因为非光滑牛顿法求解子问题要花费很多时间，其实际运算效率并不高。下面的数值例子也表明，非光滑牛顿法不一定收敛，并且算法3.1比非光滑牛顿法收敛效率更快。

Table 1. n = 100 , [ a , b ] = [ 0 , 1 ]

表1. $n=100,\left[a,b\right]=\left[0,1\right]$

Table 2. p = 10 , [ a , b ] = [ 0 , 1 ]

表2. $p=10,\left[a,b\right]=\left[0,1\right]$

Table 3. n = 100 , [ a , b ] = [ 1 , 10 ]

表3. $n=100,\left[a,b\right]=\left[1,10\right]$

Table 4. p = 10 , [ a , b ] = [ 1 , 10 ]

表4. $p=10,\left[a,b\right]=\left[1,10\right]$

在以下数值实验中，我们采用编程软件为Matlab2016a。初始矩阵 $A^0\in R^{n\times n}$ 为随机产生的n阶矩阵，iter，Err和time (s)分别表示算法的迭代次数，误差( $\text{Err}=\Vert A^{k}X-X\Lambda \Vert {}_F$ )和CPU运行时间(以秒为单位)。在每次实验时，我们随机生成一个 $n\times n$ 矩阵，其元素均匀分布在 $\left[a,b\right]$ 之间，并取其前p个特征对为给定特征对。其中第1，2个数值实验在随机生成的元素均匀分布在[0, 1]之间；第3，4个数值实验在随机生成的元素均匀分布在[1, 10]之间。具体的实验结果分别见表1~4。其中“−”表示算法迭代10,000次没有达到精度。此时，我们认为算法不收敛。

基金项目

论文得到国家自然科学基金(11661019)和贵州省数据驱动建模学习与优化创新团队(黔科合平台人才[2020]5016)资助。

NOTES

¹所有元素都是非负实数的矩阵称为非负矩阵。