求解实对称互补特征值问题的积极集方法

雷 渊，朱 琳，李 斌

（湖南大学数学学院，湖南长沙 410082）

互补特征值问题［1］（EiCP）是一类特殊的互补问题［2］，又称为锥约束下的特征值问题［3-5］，它是经典特征值问题的推广。近年来，互补特征值问题受到了广泛的关注，在结构机械系统的动力学分析、信号处理、摩擦弹性系统和声波系统的稳定性分析以及流体动力学等方面都涉及到互补特征值问题的理论和对应的特征值与特征向量的计算问题［6-8］。

目前关于互补特征值问题研究主要集中于2类常用的自对偶锥：非负锥和二阶锥［9-10］。当锥为非负锥Rn+时，该问题又称为Pareto特征值问题，其数学描述为：给定n阶实矩阵A和对称正定矩阵B，求实数λ∈R和非零向量x∈Rn满足

称满足式（1）的实数λ是矩阵对(A，B)的一个Pareto特征值，x是相应的互补特征向量。当式（1）中的矩阵A为实对称矩阵时，该问题可称为对称Pareto特征值问题。由于Pareto特征值问题（1）等价于单纯形上一个连续函数的变分不等式问题，所以该问题的解一定存在［11］。

与线性或非线性互补问题类似，利用非线性互补函数可以将Pareto特征值问题等价地转化成一个非光滑方程组，然后用半光滑牛顿方法求解［12］，但是它在许多情况下会失败。另外，求解约束规划问题的内点算法［13］、投影算法［14］以及DC（difference of convex functions）算法［15］都可用于计算大规模实对称互补特征值，但都存在计算效率不高的问题。在文献［16］中提出了一种具有谱选择搜索方向的块积极集算法（BAS），并通过大量的计算实验证明了其有效性。

序列二次规划（SQP）算法是求解非线性规划问题的一种高效算法［17-18］。序列二次规划算法通过构造一系列二次规划子问题来求解非线性规划问题，而积极集方法适合求解二次规划问题。本文主要将积极集方法和序列二次规划算法相结合，构造求解对称Pareto特征值问题的一类积极集方法。

1 求解对称Pareto特征值的序列二次规划方法

主要讨论对称Pareto特征值问题，即矩阵A、B都是对称矩阵，其中B是对称正定矩阵。当矩阵A非正定，可以找到一个足够大的正常数δ＞0，使得矩阵A+δB是正定的。若(λ，x)是矩阵对(A+δB，B)的Pareto特征值和特征向量，则由式（1）可以很容易地得到(λ−δ，x)是矩阵对(A，B)的一个Pareto特征值和特征向量。因此下面仅考虑A、B都是对称正定矩阵的情况，用符号Sn+表示n阶实对称正定矩阵的集合。

在求解对称Pareto特征值问题时，可以将原问题（1）转化为如下约束优化问题：

其拉格朗日函数为

显然问题（2）的稳定点满足如下KKT（Karush-Kuhn-Tucker）条件：

式中：λ∈R、z∈Rn+为相应等式约束和不等式约束的拉格朗日乘子。从式（4）可以知道(λ，x)是原问题（1）的一个解。因此可以利用序列二次规划方法找到优化问题（2）的一个稳定点，继而得到原问题（1）的一个Pareto特征值和其特征向量。

首先由当前迭代点x k构造二次规划子问题：

其中矩阵M是个对称正定矩阵。当可行域非空时，问题（5）有解且有唯一解，记问题（5）的解为d k。用d k作为第k次迭代的搜索方向，再选取适当步长αk，从而由迭代格式为x k+1=x k+αk d k生成序列{x k}。可以证明当问题（5）的解d k=0时，x k就是原问题（1）的一个Pareto特征向量。

下面定义一个罚函数Pσ(x)：

式中：σ是罚参数，σ＞0。在经典的序列二次规划算法中，当给定的罚参数σ对所有的k都满足σ＞|λk|，其中λk是问题（5）等式约束的拉格朗日乘子，则d k是罚函数Pσ(x)的下降方向。又因为罚函数Pσ(x)＞0有下界，故算法可收敛。在文献［19］中有类似证明，这里不再阐述。在后面的数值实验中，发现这确实可以实现，其罚参数σ随着每次迭代而更新。

简单描述用SQP算法求解对称正定矩阵Pareto特征值的基本框架如下。

（1）初始值。①给定对称正定矩阵A，B∈Sn+；②选取对称正定矩阵M∈Sn+和罚参数σ＞0；③选取误差参数ε＞0和初始迭代点x0≥0，k=0。

（2）迭代。①求出子问题（5）唯一解d k；②如果，则算法停止，输出x k。否则，转至步骤③；③选取恰当的步长αk∈(0，1]，使得Pσ(x k+αk d k)=转至步骤①。

从上面算法可知，利用序列二次规划问题求Pareto特征值问题需要找到适当的步长αk∈(0，1]，使 得在 研 究过程中发现函数Pσ(x k+αk d k)在区间(0，1]上是分段连续函数，可以直接计算出需要的步长，这里就不赘述求解过程，直接给出结果（详见文献［19］）。

定理1设x k和d k由算法SQP生成，假设d k≠0且σ＞ρ，其中ρ是矩阵B−1A的谱半径。如式（1）和式（2）选取的步长αk：

（1）当ck＞0时

（2）当ck≤0时

定理2设{x k}是算法SQP生成的迭代序列，当罚参数σ对所有的k都满足σ＞max{|λk|，ρ}，其中λk是问题（5）等式约束的拉格朗日乘子，ρ是矩阵B−1A的谱半径，则序列{x k}的任意聚点都是原问题（1）的一个Pareto特征向量。

从上面可以发现算法SQP的主要思想是用子问题（5）解d k作为搜索方向，生成迭代序列，最后收敛到原问题的解。如何快速精确求解子问题（5）是问题的关键，下面利用积极集方法求解子问题。

2 积极集算法

众所周知，积极集方法被认为是求解中小规模二次规划问题的高效的算法之一。本节主要目的是构造一类有效求解子问题（5）的解d k的积极集方法。已知当前点x k，为了简洁，将下面的量重新进行标记：，则可行域可简记为D：={d|bTd=c，d+x k≥0}。进而子问题（5）可以简写为

对于对称正定矩阵M，当可行域D非空时，易知问题（7）存在唯一解。那么d k是问题（7）的解当且仅当存在(λk，zk)满足下面KKT条件：

其中λk∈R、zk∈Rn+分别是等式约束和不等式约束的拉格朗日乘子。

下面定义2个指标集：Jk=J(d k)={i|(d k)i+(x k)i=0，i∈N}，Ik=N/Jk，称Jk为问题（7）在点d k处的积极集，其中N={1，2，…，n}。若d k是问题（7）的解，不难发现d k也是如下二次规划问题的稳定点：

引理1若问题（9）有解，则其解为

其中满足下面方程：

证明不失一般性，对矩阵和向量进行如下分块：

其中M Jk，Ik表示抽取矩阵M的Jk行和Ik列上的元素组成的主子矩阵，d Jk={d i|i∈Jk}。令

则问题（9）转化为

又因为d Jk=−x Jk，所以只需要求出的值即可。显然可以通过解下面优化问题得到要求的值。

问题（12）的解满足KKT条件：

其中是其相应的拉格朗日乘子，显然解是下面方程的解：

证毕。

综上可知，二次规划问题（7）的解也是问题（9）的解，可通过求解问题（9）来找到问题（7）的解。由于问题（9）依赖于问题（7）的最优解d k，所以解决问题（7）的关键步是确定哪些不等式约束是积极的，即设计一类确定指标集Jk的策略。

下面建立积极集法求解子问题（7）的基本框架。

（1）初始值。①给定对称正定矩阵M∈Sn+，向量x k∈Rn，q∈Rn，b∈Rn及标量c∈R；②选取初始点d0∈D，令k=0，K k为空集；③确定初始积极集：Jk=J(d 0)={i|(d 0)i+(x k)i=0}，Ik=N/Jk。

取d k+1=，K k+1为 空 集，Jk+1={i|(d k+1)i+(x k)i=0}，Ik+1=N/Jk+1；④k=k+1，转 至 步骤①。

接下来讨论积极集算法的一些基本性质。

定理3设若序列{d k}由积极集算法生成，其中d0∈D，x k≥0，则有如下性质：①对任意的k都有，d k∈D，即bTd k=c，d k+x k≥0；②对任意的k都有，g(d k+1)≤g(d k)，且等号成立当且仅当d k+1=d k。

证明（1）用归纳法证明。因为d 0∈D，不妨假设d k∈D成立，所以只需证d k+1∈D即可。在积极集算法的情况（a）、（b）中d k+1=，d k+1是问题（9）的解，显然有d k+1∈D。在情况（c）中分2种情况讨论。

情况1：当K k∩Sk=∅时，因为d k+x k≥0，所以 对∀i∈Sk，有d k+x k＞0。又 因 为Sk=所以对∀i∈Sk，有−(x k)i＜(d k)i，进而可知β*＞0。而

综上可知0＜β*＜1。令d*=d(β*)=d k+，下面证明d*∈D。

不等式约束分2种情况如下。①当i∉Sk时，易知。所以有

②当i∈Sk时，易知。所以有

综上可知d*+x k≥0，所以d*∈D。又因为d k+1=d(βk)，0≤βk≤β*，显 然d k+1也 满 足bTd k+1=c，d k+1+x k≥0，即d k+1∈D。

情况2：当K k∩Sk≠∅时，则d k+1是问题（13）的解，显然有d k+1∈D。

（2）由上可知d k、d k+1都在可行域内，由于d k和d k+1有多种取法，下面分情况讨论。

①若d k+1=，那么d k+1是问题（9）的稳定点，显然d k也满足问题（9）的约束条件，则易得g(d k+1)=即g(d k+1)≤g(d k)，等 号成立当且仅当d k+1=d k。②若d k+1=d(βk)=由 算 法 可 知g(d(βk))=显然有g(d k+1)≤g(d k)，等号成立当且仅当d k+1=d k。③若d k+1=，即d k+1是问题（13）的解，显然d k也满足问题（13）的约束条件，同①可知g(d k+1)≤g(d k)，等号成立当且仅当d k+1=d k。

当d k=且Jk=Jk−1/K k、或 者d k=d k+、或者这3种情况时，在这里就不一一讨论了，与情况且=相同，可证g(d k+1)≤g(d k)，等号成立当且仅当d k+1=d k。证毕。

从上面定理3可知由积极集算法生成的序列{d k}在可行域内，且保证目标函数g(d)下降。

定理4若是算法输出的解，则有如下性质：，综 上 可 知是问题（7）的解。

证明由积极集算法和定理3，性质①②③④显然易得，所以只需证明。因为，即

所以

定理4表明可以通过积极集算法找到问题（7）的解，即子问题（5）的解。积极集算法相当于是个组合问题，这保证了算法会在有限步结束，接下来讨论算法中一些具体问题，如在算法迭代步③中求解出βk使得

定理5在积极集算法中，若g(d(βk))=，则

证明令，则。而

所以

又因为矩阵M是正定的，则有所以g(d(β))的最小值点为。显然当βk满足式（14）时，使得证毕。

最后对问题（9）是否有解进行讨论。从引理1可以知道，当方程（10）有解时，显然问题（9）也有解。因为矩阵M是正定的，则它的主子矩阵MI也是正定矩阵，即矩阵MI可逆，那么有如下2种情况：①bIk≠0，则。系数矩阵行列式不等于零，即矩阵可逆，方程有唯一解。②bIk=0，则方程无解或有无穷多个解。

为了避免第2种情况出现，即保证bIk≠0，那么可对指标集Jk做如下处理。因为b为非零向量，则必可找到一个i使得bi≠0。若i∈Jk，则Jk=Jk/i；若i∉Jk，则指标集Jk保持不变。这样保证了b Ik≠0，方程（10）有唯一解，所以在实际计算中可以保证问题（9）有解。

3 数值实验

本节详细描述数值算例的设计，尝试通过数值实验验证积极集方法的有效性以及其在计算时间、迭代步数以及互补精度方面均优于已有算法。所有数值实验在操作系统为64位windows 10和处理器为Intel（R）Core（TM）i7-8700 CPU@3.20GHz 3.19 GHz的电脑上进行，使用的软件是Matlab R2018b。

SQP算法最不可缺少的部分是子问题（5）中对称正定矩阵M的选取。因为本文主要讨论矩阵A是对称正定的情况，在这里直接令M=A。为了简洁，用积极集算法求解子问题的序列二次规划问题，记为SQP（J）。又因为Matlab软件中的内置函数“quadprog”能直接求解二次规划，在所有实验中算法SQP（Q）表示直接用内置函数“quadprog”求解子问题（5）。

积极集算法中求解子问题的初始向量d0的公式如下：

其中，i为Bx0中第1个分量为正的指标。因为矩阵B为正定矩阵，这样的分量一定存在。这样选取保证了d0在可行域内，也从侧面证明了子问题可行域非空。在第k次外迭代后，积极集算法中求解子问题的初始向量为d k=(1−αk−1)d k−1，其中αk−1和d k−1分别为上一步外迭代得到的步长和子问题的解。另外，积极集算法在迭代步③中（b）的情况2中m的选取也会影响迭代次数，在这里取，其中l表示小于零的分量的个数。

为了有效地分析算法SQP（J）和SQP（Q）性能，也用算法BAS［16］求解了所有数值实验。在本节的所有数值实验中，初始矩阵A、B都为随机生成的对称正定矩阵。初始向量为x 0=es，其中es表示第s个分量为1、其余分量为零的向量，且s=argmax{r i|i=1，…，n}，r i=min{ajibii−aiibji|j=1，…，n}。所 有算法的停止准则为。另外，还讨论了算法的互补性，即

例1为了测试算法的可行性，首先计算一个简单的问题。在这里令

初始向量为x0=[1 0 1 0]T，通过算法SQP（J）和算法SQP（Q）得到

这意味这λ*和x*是矩阵对(A，B)的Pareto特征值和对应的Pareto特征向量。

例2为了进一步测试算法的有效性，测试一些随机生成的问题。随机生成5组阶数分别为50、500和1 000的对称正定矩阵对(A，B)，分别列出了算法BAS、算法SQP（Q）和算法SQP（J）计算矩阵对(A，B)的Pareto特征值所需的CPU运行时间（s），迭代次数IT以及互补性H，见表1-3。其中，N表示矩阵的阶数。

从表1、表2和表3可以看出在阶数N=50时，算法BAS的运行时间少于算法SQP（J），但互补性精度低于算法SQP（J）。当矩阵阶数增高到N=1000时，算法SQP（J）的运行时间明显低于算法BAS，而算法SQP（Q）的互补性表现糟糕。从整体运算时间和互补性上来看，算法SQP（J）比算法BAS和算法SQP（Q）性能更稳，效率更高。

表1 随机生成的50阶问题的数值结果比较Tab.1 Comparison of numerical results of randomly generated 50th order problems

表2 随机生成的500阶问题的数值结果比较Tab.2 Comparison of numerical results of randomly generated 500th order problems

表3 随机生成的1 000阶问题的数值结果比较Tab.3 Comparison of numerical results of randomly generated 1000th order problems

4 结论

构造了求解实对称互补特征值问题的积极集方法。实验数值结果表明，所提出的算法SQP（J）在迭代步数和计算时间以及互补性上均优于算法BAS和算法SQP（Q）。算法SQP（J）能够快速求解实对称互补特征值问题，主要原因在于积极集算法更加适合求解子问题。如何选取子问题中的矩阵M是接下来需要研究的问题。

作者贡献声明：

雷 渊：提出研究思路，设计研究方案，修订最终版本。

朱 琳：完成研究方案，进行数值实验，分析数据，起草论文。

李 斌：负责数值实验的编程。