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

摘要

基于序列二次规划算法构造了求解实对称互补特征值问题的一类积极集方法。 通过特殊的积极集指标选取策略，该积极集方法计算得到的迭代序列具有单调下降特征，并从理论上证明了该方法的收敛性。 数值实验结果表明该方法是行之有效的，并且在互补性和迭代时间上均优于Matlab软件的内置算法。

关键词

互补特征值问题; 非负锥; 序列二次规划; 积极集

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

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

$$\mathit{x}\ge \mathrm{0},\mathit{A}\mathit{x}-\lambda \mathit{B}\mathit{x}\ge \mathrm{0},{\mathit{x}}^{\mathrm{T}}(\mathit{A}\mathit{x}-\lambda \mathit{B}\mathit{x})=\mathrm{0}$$

（1）

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

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

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

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

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

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

$$\begin{array}{c}\underset{\mathit{x}\in {\mathbf{R}}^{\mathit{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}f(\mathit{x})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{T}}\mathit{A}\mathit{x}\\ \begin{array}{c}\mathrm{s}.\mathrm{t}.\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{T}}\mathit{B}\mathit{x}-\mathrm{1}=\mathrm{0}\\ \mathit{x}\ge \mathrm{0}\end{array}\end{array}$$

（2）

其拉格朗日函数为

$$L(\mathit{x},\lambda ,\mathit{z})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{T}}\mathit{A}\mathit{x}-\lambda \left(\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{T}}\mathit{B}\mathit{x}-\mathrm{1}\right)-{\mathit{z}}^{\mathrm{T}}\mathit{x}$$

（3）

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

$$\left\{\begin{array}{l}\mathit{A}\mathit{x}-\lambda \mathit{B}\mathit{x}=\mathit{z}\\ \frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{T}}\mathit{B}\mathit{x}-\mathrm{1}=\mathrm{0}\\ {\mathit{x}}^{\mathrm{T}}\mathit{z}=\mathrm{0}\\ \mathit{z},\mathit{x}\ge \mathrm{0}\end{array}\right.$$

（4）

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

首先由当前迭代点${\mathit{x}}_{\mathit{k}}$构造二次规划子问题：

$$\begin{array}{c}\underset{\mathit{d}\in {\mathbf{R}}^{\mathit{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\mathit{M}\mathit{d}+(\mathit{A}{\mathit{x}}_{\mathit{k}}{)}^{\mathrm{T}}\mathit{d}\\ \mathrm{s}.\mathrm{t}.(\mathit{B}{\mathit{x}}_{\mathit{k}}{)}^{\mathrm{T}}\mathit{d}+\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}_{\mathit{k}}^{\mathrm{T}}\mathit{B}{\mathit{x}}_{\mathit{k}}-\mathrm{1}=\mathrm{0}\\ \mathit{d}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}\end{array}$$

（5）

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

下面定义一个罚函数$P_{\sigma }(\mathit{x})$：

$$P_{\sigma }(\mathit{x})=f(\mathit{x})+\sigma \left|\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}^{\mathrm{T}}\mathit{B}\mathit{x}-\mathrm{1}\right|$$

（6）

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

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

（1）初始值。①给定对称正定矩阵$\mathit{A},\mathit{B}\in {\mathit{S}}_{+}^{\mathit{n}}$；②选取对称正定矩阵$\mathit{M}\in {\mathit{S}}_{+}^{\mathit{n}}$和罚参数$\sigma >\mathrm{0}$；③选取误差参数$\epsilon >\mathrm{0}$和初始迭代点${\mathit{x}}_{\mathrm{0}}\ge \mathrm{0}$，$k=\mathrm{0}$。

（2）迭代。①求出子问题（5）唯一解${\mathit{d}}_{\mathit{k}}$；②如果$‖{\mathit{d}}_{\mathit{k}}‖\le \epsilon$，则算法停止，输出${\mathit{x}}_{\mathit{k}}$。否则，转至步骤③；③选取恰当的步长$\alpha {}_k{}^{}\in (\mathrm{0,1}]$，使得$P_{\sigma }({\mathit{x}}_{\mathit{k}}+{\alpha }_k{\mathit{d}}_{\mathit{k}})=\underset{\mathrm{0}<\alpha \le \mathrm{1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{P}_{\sigma }({\mathit{x}}_{\mathit{k}}+\alpha {\mathit{d}}_{\mathit{k}})$；④${\mathit{x}}_{\mathit{k}+\mathrm{1}}={\mathit{x}}_{\mathit{k}}+{\alpha }_k{\mathit{d}}_{\mathit{k}}$，$k=k+\mathrm{1}$，转至步骤①。

从上面算法可知，利用序列二次规划问题求Pareto特征值问题需要找到适当的步长$\alpha {}_k{}^{}\in (\mathrm{0,1}]$，使得$P_{\sigma }({\mathit{x}}_{\mathit{k}}+{\alpha }_k{\mathit{d}}_{\mathit{k}})=\underset{\mathrm{0}<\alpha \le \mathrm{1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{P}_{\sigma }({\mathit{x}}_{\mathit{k}}+\alpha {\mathit{d}}_{\mathit{k}})$。在研究过程中发现函数$P_{\sigma }({\mathit{x}}_{\mathit{k}}+{\alpha }_k{\mathit{d}}_{\mathit{k}})$在区间$(\mathrm{0,1}]$上是分段连续函数，可以直接计算出需要的步长，这里就不赘述求解过程，直接给出结果（详见文献［19］）。

定理 1  　设${\mathit{x}}_{\mathit{k}}$和${\mathit{d}}_{\mathit{k}}$由算法SQP生成，假设${\mathit{d}}_{\mathit{k}}\ne \mathrm{0}$且$\sigma >\rho$，其中$\rho$是矩阵${\mathit{B}}^{-\mathrm{1}}\mathit{A}$的谱半径。如式（1）和式（2）选取的步长${\alpha }_k$：

（1）当$c_k>\mathrm{0}$时

$${\alpha }_k=\left\{\begin{array}{c}\mathrm{1},\mathrm{1}\le {\delta }_k\\ {\delta }_k,{\delta }_k^{\mathrm{1}}<{\delta }_k<\mathrm{1}\\ {\delta }_k^{\mathrm{1}},{\delta }_k\le {\delta }_k^{\mathrm{1}}\end{array}\right.$$

（2）当$c_k\le \mathrm{0}$时

$${\alpha }_k=\left\{\begin{array}{c}\mathrm{1},\mathrm{1}\le {\delta }_k\\ {\delta }_k,{\delta }_k<\mathrm{1}\end{array}\right.$$

则有$P_{\sigma }({\mathit{x}}_{\mathit{k}}+{\alpha }_k{\mathit{d}}_{\mathit{k}})=\underset{\mathrm{0}<\alpha \le \mathrm{1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{P}_{\sigma }({\mathit{x}}_{\mathit{k}}+{\alpha }_{}{\mathit{d}}_{\mathit{k}})$。式中：$c_k=\mathrm{1}-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}_{\mathit{k}}^{\mathrm{T}}\mathit{B}{\mathit{x}}_{\mathit{k}}$，${\delta }_k=-\frac{{\mathit{x}}_{\mathit{k}}^{\mathrm{T}}\mathit{A}{\mathit{d}}_{\mathit{k}}+\sigma c_k}{{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}(\mathit{A}+\sigma \mathit{B}){\mathit{d}}_{\mathit{k}}}$，${\delta }_k^{\mathrm{1}}=\frac{-c_k+\sqrt[]{c_k^{\mathrm{2}}+\mathrm{2}{c}_k{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{B}{\mathit{d}}_{\mathit{k}}}}{{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{B}{\mathit{d}}_{\mathit{k}}}$。

定理 2  　设$\{x_k\}$是算法SQP生成的迭代序列，当罚参数$\sigma$对所有的$k$都满足$\sigma >\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left.|{\lambda }_k|,\rho \right\}\right.$，其中${\lambda }_k$是问题（5）等式约束的拉格朗日乘子，$\rho$是矩阵${\mathit{B}}^{-\mathrm{1}}\mathit{A}$的谱半径，则序列$\{{\mathit{x}}_{\mathit{k}}\}$的任意聚点都是原问题（1）的一个Pareto特征向量。

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

2 积极集算法

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

$$\begin{array}{c}\underset{\mathit{d}\in {\mathbf{R}}^{\mathit{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}g(\mathit{d})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\mathbf{M}\mathbf{d}+{\mathit{q}}^{\mathrm{T}}\mathit{d}\\ \mathrm{s}.\mathrm{t}.{\mathit{b}}^{\mathrm{T}}\mathit{d}=c\\ \mathit{d}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}\end{array}$$

（7）

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

$$\left\{\begin{array}{l}\mathit{M}{\mathit{d}}_{\mathit{k}}+\mathit{q}={\lambda }_k\mathit{b}+{\mathit{z}}_{\mathit{k}}\\ {\mathit{b}}^{\mathrm{T}}{\mathit{d}}_{\mathit{k}}=c\\ {\mathit{z}}_{\mathit{k}}^{\mathrm{T}}({\mathit{d}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}})=\mathrm{0}\\ {\mathit{d}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}\\ {\mathit{z}}_{\mathit{k}}\ge \mathrm{0}\end{array}\right.$$

（8）

其中${\lambda }_k\in \mathbf{R}$、${\mathit{z}}_{\mathit{k}}\in {\mathit{R}}_{+}^{\mathit{n}}$分别是等式约束和不等式约束的拉格朗日乘子。

下面定义2个指标集：$J_k=J({\mathit{d}}_{\mathit{k}})=\{i|({\mathit{d}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=\mathrm{0},i\in \mathbf{N}\},I_k=\mathbf{N}/J_k,$称$J_k$为问题（7）在点${\mathit{d}}_{\mathit{k}}$处的积极集，其中$\mathbf{N}=\{\mathrm{1,2},\cdots ,n\}$。若${\mathit{d}}_{\mathit{k}}$是问题（7）的解，不难发现${\mathit{d}}_{\mathit{k}}$也是如下二次规划问题的稳定点：

$$\begin{array}{c}\underset{\mathit{d}\in {\mathbf{R}}^{\mathit{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\widehat{g}(\mathit{d})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\mathit{M}\mathit{d}+{\mathit{q}}^{\mathrm{T}}\mathit{d}\\ \mathrm{s}.\mathrm{t}.{\mathit{b}}^{\mathrm{T}}\mathit{d}=c\\ d_i=-({\mathit{x}}_{\mathit{k}}{)}_i,i\in J_k\end{array}$$

（9）

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

$${\widehat{\mathit{d}}}_{\mathit{k}}=\left(\begin{array}{c}-{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}\\ {\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)$$

其中${\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}$满足下面方程：

$$\left[\begin{array}{cc}{\mathit{M}}_{\mathit{I}}& {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}^{\mathrm{T}}& \mathrm{0}\end{array}\right]\left(\begin{array}{c}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}\\ -{\widehat{\lambda }}_k\end{array}\right)=\left(\begin{array}{c}{\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}-{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}+c\end{array}\right)$$

（10）

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

$$\mathit{M}=\left[\begin{array}{cc}{\mathit{M}}_{{\mathit{J}}_{\mathit{k}},{\mathit{J}}_{\mathit{k}}}& {\mathit{M}}_{{\mathit{J}}_{\mathit{k}},{\mathit{I}}_{\mathit{k}}}\\ {\mathit{M}}_{{\mathit{I}}_{\mathit{k}},{\mathit{J}}_{\mathit{k}}}& {\mathit{M}}_{{\mathit{I}}_{\mathit{k}},{\mathit{I}}_{\mathit{k}}}\end{array}\right]=\left[\begin{array}{cc}{\mathit{M}}_{\mathit{J}}& {\mathit{M}}_{\mathit{J},\mathit{I}}\\ {\mathit{M}}_{\mathit{I},\mathit{J}}& {\mathit{M}}_{\mathit{I}}\end{array}\right]$$

$$\mathit{b}=\left(\begin{array}{c}{\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right),\mathit{q}=\left(\begin{array}{c}{\mathit{q}}_{{\mathit{J}}_{\mathit{k}}}\\ {\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)$$

$$\mathit{d}=\left(\begin{array}{c}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}\\ {\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right),{\mathit{x}}_{\mathit{k}}=\left(\begin{array}{c}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}\\ {\mathit{x}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)$$

其中${\mathit{M}}_{{\mathit{J}}_{\mathit{k}},{\mathit{I}}_{\mathit{k}}}$表示抽取矩阵$\mathit{M}$的$J_k$行和$I_k$列上的元素组成的主子矩阵，${\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}=\{{\mathit{d}}_{\mathit{i}}|i\in J_k\}$。令

$$\begin{array}{c}\widehat{g}(\mathit{d})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\mathbf{M}\mathbf{d}+{\mathit{q}}^{\mathrm{T}}\mathit{d}=\\ \begin{array}{l}\frac{\mathrm{1}}{\mathrm{2}}{{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{M}}_{\mathit{I}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}+({\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}+{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}{)}^{\mathrm{T}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}+\\ \frac{\mathrm{1}}{\mathrm{2}}{{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{M}}_{\mathit{J}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}+{{\mathit{q}}_{{\mathit{J}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}\end{array}\end{array}$$

则问题（9）转化为

$$\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}\widehat{g}(\mathit{d})=\frac{\mathrm{1}}{\mathrm{2}}{{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{M}}_{\mathit{I}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}+({\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}+\\ {\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}{)}^{\mathrm{T}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}+\frac{\mathrm{1}}{\mathrm{2}}{{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{M}}_{\mathit{J}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}+{{\mathit{q}}_{{\mathit{J}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}\\ \mathrm{s}.\mathrm{t}.{\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}=-{\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}+c\\ {\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}=-{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}\end{array}$$

（11）

又因为${\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}=-{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}$，所以只需要求出${\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}$的值即可。显然可以通过解下面优化问题得到要求的值${\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}$。

$$\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}\overline{g}({\mathit{d}}_{{\mathit{I}}_{\mathit{k}}})=\frac{\mathrm{1}}{\mathrm{2}}{{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}}^{\mathrm{T}}{\mathit{M}}_{\mathit{I}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}+(-{\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}+{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}{)}^{\mathrm{T}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}\\ \mathrm{s}.\mathrm{t}.{\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}={\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}+c\end{array}$$

（12）

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

$$\left\{\begin{array}{l}{\mathit{M}}_{\mathit{I}}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}-{\widehat{\lambda }}_{\mathit{k}}{\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}={\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}-{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}^{\mathrm{T}}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}={\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}+c\end{array}\right.$$

其中${\widehat{\lambda }}_k$是其相应的拉格朗日乘子，显然解$({\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}},{\widehat{\lambda }}_k)$是下面方程的解：

$$\left[\begin{array}{cc}{\mathit{M}}_{\mathit{I}}& {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}^{\mathrm{T}}& \mathrm{0}\end{array}\right]\left(\begin{array}{c}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}\\ -{\widehat{\lambda }}_k\end{array}\right)=\left(\begin{array}{c}{\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}-{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{J}}_{\mathit{K}}}^{\mathrm{T}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}+c\end{array}\right)$$

证毕。

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

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

（1）初始值。①给定对称正定矩阵$\mathit{M}\in$${\mathit{S}}_{+}^{\mathit{n}}$，向量${\mathit{x}}_{\mathit{k}}\in {\mathit{R}}^{\mathit{n}}$，$\mathit{q}\in {\mathit{R}}^{\mathit{n}}$，$\mathit{b}\in {\mathit{R}}^{\mathit{n}}$及标量$c\in \mathbf{R}$；②选取初始点${\mathit{d}}_{\mathrm{0}}\in \mathit{D}$，令$k=\mathrm{0}$，$K_k$为空集；③确定初始积极集：$J_k=J({\mathit{d}}_{\mathrm{0}})=\{i|({\mathit{d}}_{\mathrm{0}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=\mathrm{0}\}$，$I_k=\mathbf{N}/J_k$。

（2）迭代。①求问题（9）的解$({\widehat{\mathit{d}}}_{\mathit{k}},{\widehat{\lambda }}_k)$；②求${\widehat{\mathit{z}}}_{\mathit{k}}=\mathit{M}{\widehat{\mathit{d}}}_{\mathit{k}}+\mathit{q}-{\widehat{\lambda }}_k\mathit{b}$；③判断。（a）若${\widehat{\mathit{d}}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$，且${\widehat{\mathit{z}}}_{\mathit{k}}\ge \mathrm{0}$，则算法停止，输出$({\widehat{\mathit{d}}}_{\mathit{k}},{\widehat{\lambda }}_k,{\widehat{\mathit{z}}}_{\mathit{k}})$；（b）若${\widehat{\mathit{d}}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$，但存在指标$i$使得$({\widehat{\mathit{z}}}_{\mathit{k}}{)}_i<\mathrm{0}$，则令${\widehat{J}}_k=\{i|({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=\mathrm{0}\}$，显然有$J_k\subseteq {\widehat{J}}_k$，这里又分2种情况讨论：情况1，若$J_k\subset {\widehat{J}}_k$，则取${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\widehat{\mathit{d}}}_{\mathit{k}}$，$K_{k+\mathrm{1}}$为空集，$J_{k+\mathrm{1}}={\widehat{J}}_k$，$I_{k+\mathrm{1}}=\mathbf{N}/J_{k+\mathrm{1}}$；情况2，若$J_k={\widehat{J}}_k$，取${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\widehat{\mathit{d}}}_{\mathit{k}}$，$K_{k+\mathrm{1}}=\{i\in {\widehat{J}}_k|({\widehat{\mathit{z}}}_{\mathit{k}}{)}_i<({\widehat{\mathit{z}}}_{\mathit{k}}{)}_m\}$，$J_{k+\mathrm{1}}=J_k/K_{k+\mathrm{1}}$，$I_{k+\mathrm{1}}=N/J_{k+\mathrm{1}}$，其中${\left({\widehat{\mathit{z}}}_{\mathit{k}}\right)}_m$表示向量${\widehat{\mathit{z}}}_{\mathit{k}}$中第$m$小的分量；（c）若存在$i$使$({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i<\mathrm{0}$，即${\widehat{\mathit{d}}}_{\mathit{k}}\notin \mathit{D}$。令$S_k=\{i\in I_k|({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i<\mathrm{0}\}$，分2种情况讨论：情况1，若$K_k\bigcap S_k=\mathrm{\varnothing }$，令${\beta }^{\mathrm{*}}=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\frac{({\mathit{x}}_{\mathit{k}}{)}_i+({\mathit{d}}_{\mathit{k}}{)}_i}{({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}|i\in S_k\right\}$，$\mathit{d}(\beta )={\mathit{d}}_{\mathit{k}}+\beta ({\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}})$，求出${\beta }_k$使得$g(\mathit{d}({\beta }_k))=\underset{\mathrm{0}\le {\beta }_k\le {\beta }^{\mathrm{*}}}{\mathrm{m}\mathrm{i}\mathrm{n}}g(\mathit{d}(\beta ))$，则取${\mathit{d}}_{\mathit{k}+\mathrm{1}}=\mathit{d}\left({\beta }_k\right)$，$K_{k+\mathrm{1}}$为空集，$J_{k+\mathrm{1}}=\{i|({\mathit{d}}_{\mathit{k}+\mathrm{1}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=\mathrm{0}\}$，$I_{k+\mathrm{1}}=\mathbf{N}/J_{k+\mathrm{1}}$；情况2，若$K_k\bigcap S_k\ne \mathrm{\varnothing }$，求出下面优化问题的解${\tilde{\mathit{d}}}_{\mathit{k}}$：

$$\begin{array}{c}\underset{\mathrm{d}\in {\mathbf{R}}^{\mathrm{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\tilde{g}(\mathit{d})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\mathbf{M}\mathbf{d}+{\mathit{q}}^{\mathrm{T}}\mathit{d}\\ \begin{array}{l}\mathrm{s}.\mathrm{t}.{\mathit{b}}^{\mathrm{T}}\mathit{d}=c\\ {\mathit{d}}_{{\mathit{J}}_{\mathit{k}}}=-({\mathit{x}}_{\mathit{k}}{)}_{{\mathit{J}}_{\mathit{k}}}\end{array}\\ {\mathit{d}}_{{\mathit{I}}_{\mathit{k}}}+({\mathit{x}}_{\mathit{k}}{)}_{{\mathit{I}}_{\mathit{k}}}\ge \mathrm{0}\end{array}$$

（13）

取${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\tilde{\mathit{d}}}_{\mathit{k}}$，$K_{k+\mathrm{1}}$为空集，$J_{k+\mathrm{1}}=\{i|({\mathit{d}}_{\mathit{k}+\mathrm{1}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=\mathrm{0}\}$，$I_{k+\mathrm{1}}=\mathbf{N}/J_{k+\mathrm{1}}$；④$k=k+\mathrm{1}$，转至步骤①。

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

定理 3  　设若序列$\{{\mathit{d}}_{\mathit{k}}\}$由积极集算法生成，其中${\mathit{d}}_{\mathrm{0}}\in \mathit{D}$，${\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$，则有如下性质：①对任意的$k$都有，${\mathit{d}}_{\mathit{k}}\in \mathit{D}$，即${\mathit{b}}^{\mathrm{T}}{\mathit{d}}_{\mathit{k}}=c$，${\mathit{d}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$；②对任意的$k$都有，$g({\mathit{d}}_{\mathit{k}+\mathrm{1}})\le g({\mathit{d}}_{\mathit{k}})$，且等号成立当且仅当${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\mathit{d}}_{\mathit{k}}$。

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

情况1：当$K_k\bigcap S_k=\mathrm{\varnothing }$时，因为${\mathit{d}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$，所以对$\forall i\in S_k$，有${\mathit{d}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}>\mathrm{0}$。又因为$S_k=\{i\in I_k|({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i<\mathrm{0}\}$，所以对$\forall i\in S_k$，有$({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i<-{\left({\mathit{x}}_{\mathit{k}}\right)}_i<{\left({\mathit{d}}_{\mathit{k}}\right)}_i$，进而可知${\beta }^{\mathrm{*}}>\mathrm{0}$。而

$$\begin{array}{l}\frac{({\mathit{x}}_{\mathit{k}}{)}_i+({\mathit{d}}_{\mathit{k}}{)}_i}{({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}=\frac{({\mathit{x}}_{\mathit{k}}{)}_i+({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}{({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}=\\ \mathrm{1}+\frac{({\mathit{x}}_{\mathit{k}}{)}_i+({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}{({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}<\mathrm{1}\end{array}$$

综上可知$\mathrm{0}<{\beta }^{\mathrm{*}}<\mathrm{1}$。令${\mathit{d}}^{\mathit{*}}=\mathit{d}({\beta }^{\mathrm{*}})={\mathit{d}}_{\mathit{k}}+{\beta }^{\mathrm{*}}({\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}})$，下面证明${\mathit{d}}^{\mathit{*}}\in \mathit{D}$。

等式约束：${\mathit{b}}^{\mathrm{T}}{\mathit{d}}^{\mathit{*}}={\mathit{b}}^{\mathrm{T}}({\mathit{d}}_{\mathit{k}}+{\beta }^{\mathrm{*}}({\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}}))={\mathit{b}}^{\mathrm{T}}{\mathit{d}}_{\mathit{k}}+{\beta }^{\mathrm{*}}({\mathit{b}}^{\mathrm{T}}{\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{b}}^{\mathrm{T}}{\mathit{d}}_{\mathit{k}})=c$。

不等式约束分2种情况如下。①当$i\notin S_k$时，易知$({\mathit{d}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i\ge \mathrm{0}$，$({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i\ge \mathrm{0}$。所以有

$$\begin{array}{c}\begin{array}{l}({\mathit{d}}^{\mathit{*}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=({\mathit{d}}_{\mathit{k}}{)}_i+{\beta }^{\mathrm{*}}(({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i-({\mathit{d}}_{\mathit{k}}{)}_i)+\\ ({\mathit{x}}_{\mathit{k}}{)}_i=(\mathrm{1}-{\beta }^{\mathrm{*}})({\mathit{d}}_{\mathit{k}}{)}_i+\end{array}\\ {\beta }^{\mathrm{*}}({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}})i=(\mathrm{1}-{\beta }^{\mathrm{*}})(({\mathit{d}}_{\mathit{k}}{)}_i{+}_{}\\ ({\mathit{x}}_{\mathit{k}}{)}_i)+{\beta }^{\mathrm{*}}(({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i)\ge \mathrm{0}\\ \end{array}$$

②当$i\in S_k$时，易知$({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i-{\left({\mathit{d}}_{\mathit{k}}\right)}_i<\mathrm{0}$。所以有

$$\begin{array}{c}({\mathit{d}}^{\mathit{*}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=({\mathit{d}}_{\mathit{k}}{)}_i+{\beta }^{\mathrm{*}}(({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i-({\mathit{d}}_{\mathit{k}}{)}_i)+\\ ({\mathit{x}}_{\mathit{k}}{)}_i=({\mathit{d}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i+{\beta }^{\mathrm{*}}(({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i-({\mathit{d}}_{\mathit{k}}{)}_i)\ge \\ ({\mathit{d}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i+\frac{({\mathit{d}}_{\mathit{k}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i}{({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}(({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i\\ -({\mathit{d}}_{\mathit{k}}{)}_i)=\mathrm{0}\end{array}$$

综上可知${\mathit{d}}^{\mathit{*}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$，所以${\mathit{d}}^{\mathit{*}}\in \mathit{D}$。又因为${\mathit{d}}_{\mathit{k}+\mathrm{1}}=\mathit{d}\left({\beta }_k\right)$，$\mathrm{0}\le {\beta }_k\le {\beta }^{\mathrm{*}}$，显然${\mathit{d}}_{\mathit{k}+\mathrm{1}}$也满足${\mathit{b}}^{\mathrm{T}}{\mathit{d}}_{\mathit{k}+\mathrm{1}}=c$，${\mathit{d}}_{\mathit{k}+\mathrm{1}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$，即${\mathit{d}}_{\mathit{k}+\mathrm{1}}\in \mathit{D}$。

情况2：当$K_k\bigcap S_k\ne \mathrm{\varnothing }$时，则${\mathit{d}}_{\mathit{k}+\mathrm{1}}$是问题（13）的解，显然有${\mathit{d}}_{\mathit{k}+\mathrm{1}}\in \mathit{D}$。

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

当${\mathit{d}}_{\mathit{k}}={\widehat{\mathit{d}}}_{\mathit{k}-\mathrm{1}}$时，且${\widehat{J}}_k=\{i|({\widehat{\mathit{d}}}_{\mathit{k}-\mathrm{1}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=$$\mathrm{0}\}$，即$d_k$是如下问题的解：

$$\begin{array}{l}\underset{\mathit{d}\in {\mathbf{R}}^{\mathit{n}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\widehat{g}(\mathit{d})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}^{\mathrm{T}}\mathit{M}\mathit{d}+{\mathit{q}}^{\mathrm{T}}\mathit{d}\\ \mathrm{s}.\mathrm{t}.{\mathit{b}}^{\mathrm{T}}\mathit{d}=c\\ {\mathit{d}}_i=-({\mathit{x}}_{\mathit{k}}{)}_i\\ i\in J_{k-\mathrm{1}}\end{array}$$

①若${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\widehat{\mathit{d}}}_{\mathit{k}}$，那么${\mathit{d}}_{\mathit{k}+\mathrm{1}}$是问题（9）的稳定点，显然${\mathit{d}}_{\mathit{k}}$也满足问题（9）的约束条件，则易得$g({\mathit{d}}_{\mathit{k}+\mathrm{1}})=\widehat{g}({\mathit{d}}_{\mathit{k}+\mathrm{1}})\le \widehat{g}({\mathit{d}}_{\mathit{k}})=g({\mathit{d}}_{\mathit{k}})$，即$g({\mathit{d}}_{\mathit{k}+\mathrm{1}})\le g({\mathit{d}}_{\mathit{k}})$，等号成立当且仅当${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\mathit{d}}_{\mathit{k}}$。②若${\mathit{d}}_{\mathit{k}+\mathrm{1}}=\mathit{d}({\beta }_k)={\mathit{d}}_{\mathit{k}}+{\beta }_k({\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}})$，由算法可知$g(\mathit{d}({\beta }_k))=\underset{\mathrm{0}\le {\beta }_k\le {\beta }^{\mathrm{*}}}{\mathrm{m}\mathrm{i}\mathrm{n}}g(\mathit{d}(\beta ))$，显然有$g({\mathit{d}}_{\mathit{k}+\mathrm{1}})\le g({\mathit{d}}_{\mathit{k}})$，等号成立当且仅当${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\mathit{d}}_{\mathit{k}}$。③若${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\tilde{\mathit{d}}}_{\mathit{k}}$，即${\mathit{d}}_{\mathit{k}+\mathrm{1}}$是问题（13）的解，显然${\mathit{d}}_{\mathit{k}}$也满足问题（13）的约束条件，同①可知$g({\mathit{d}}_{\mathit{k}+\mathrm{1}})\le g({\mathit{d}}_{\mathit{k}})$，等号成立当且仅当${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\mathit{d}}_{\mathit{k}}$。

当${\mathit{d}}_{\mathit{k}}={\widehat{\mathit{d}}}_{\mathit{k}-\mathrm{1}}$且$J_k=J_{k-\mathrm{1}}/K_k$、或者${\mathit{d}}_{\mathit{k}}={\mathit{d}}_{\mathit{k}}+{\beta }_{k-\mathrm{1}}({\widehat{\mathit{d}}}_{\mathit{k}-\mathrm{1}}-{\mathit{d}}_{\mathit{k}-\mathrm{1}})$、或者${\mathit{d}}_{\mathit{k}}={\tilde{\mathit{d}}}_{\mathit{k}-\mathrm{1}}$这3种情况时，在这里就不一一讨论了，与情况${\mathit{d}}_{\mathit{k}}={\widehat{\mathit{d}}}_{\mathit{k}-\mathrm{1}}$且${\widehat{J}}_k=\{i|({\widehat{\mathit{d}}}_{\mathit{k}-\mathrm{1}}{)}_i+({\mathit{x}}_{\mathit{k}}{)}_i=\mathrm{0}\}$相同，可证$g({\mathit{d}}_{\mathit{k}+\mathrm{1}})\le g({\mathit{d}}_{\mathit{k}})$，等号成立当且仅当${\mathit{d}}_{\mathit{k}+\mathrm{1}}={\mathit{d}}_{\mathit{k}}$。证毕。

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

定理 4  　若$({\widehat{\mathit{d}}}_{\mathit{k}},{\widehat{\lambda }}_k,{\widehat{\mathit{z}}}_{\mathit{k}})$是算法输出的解，则有如下性质：①${\mathit{b}}^{\mathrm{T}}{\widehat{\mathit{d}}}_{\mathit{k}}=c$；②${\widehat{\mathit{d}}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}}\ge \mathrm{0}$；③${\widehat{\mathit{z}}}_{\mathit{k}}=\mathit{M}{\widehat{\mathit{d}}}_{\mathit{k}}+\mathit{q}-{\widehat{\lambda }}_k\mathit{b}$；④${\widehat{\mathit{z}}}_{\mathit{k}}\ge \mathrm{0}$；⑤${\widehat{\mathit{z}}}_{\mathit{k}}^{\mathrm{T}}({\widehat{\mathit{d}}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}})=\mathrm{0}$，综上可知$({\widehat{\mathit{d}}}_{\mathit{k}},{\widehat{\lambda }}_k,{\widehat{\mathit{z}}}_{\mathit{k}})$是问题（7）的解。

证明  　由积极集算法和定理3，性质①②③④显然易得，所以只需证明${\widehat{\mathit{z}}}_{\mathit{k}}^{\mathrm{T}}({\widehat{\mathit{d}}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}})=\mathrm{0}$。因为${\widehat{\mathit{z}}}_{\mathit{k}}=\mathit{M}{\widehat{\mathit{d}}}_{\mathit{k}}+\mathit{q}-{\widehat{\lambda }}_k\mathit{b}$，即

$$\left(\begin{array}{c}{\widehat{\mathit{z}}}_{{\mathit{J}}_{\mathit{k}}}\\ {\widehat{\mathit{z}}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)=\left(\begin{array}{cc}{\mathit{M}}_{\mathit{J}}& {\mathit{M}}_{\mathit{J},\mathit{I}}\\ {\mathit{M}}_{\mathit{I},\mathit{J}}& {\mathit{M}}_{\mathit{I}}\end{array}\right)\left(\begin{array}{c}{\widehat{\mathit{d}}}_{{\mathit{J}}_{\mathit{k}}}\\ {\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)+\left(\begin{array}{c}{\mathit{q}}_{{\mathit{J}}_{\mathit{k}}}\\ {\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)-{\widehat{\lambda }}_k\left(\begin{array}{c}{\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)$$

所以

$${\widehat{\mathit{z}}}_{{\mathit{J}}_{\mathit{k}}}={\mathit{M}}_{\mathit{J}}{\widehat{\mathit{d}}}_{{\mathit{J}}_{\mathit{k}}}+{\mathit{M}}_{\mathit{J},\mathit{I}}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}+{\mathit{q}}_{{\mathit{J}}_{\mathit{k}}}-{\widehat{\lambda }}_k{\mathit{b}}_{{\mathit{J}}_{\mathit{k}}}$$

$${\widehat{\mathit{z}}}_{{\mathit{I}}_{\mathit{k}}}={\mathit{M}}_{\mathit{I},\mathit{J}}{\widehat{\mathit{d}}}_{{\mathit{J}}_{\mathit{k}}}+{\mathit{M}}_{\mathit{I}}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}+{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}-{\widehat{\lambda }}_k{\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}$$

由引理1知，${\widehat{\mathit{d}}}_{\mathit{k}}=\left(\begin{array}{c}-{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}\\ {\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}\end{array}\right)$，且$({\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}},{\widehat{\lambda }}_k)$满足KKT条件${\mathit{M}}_{\mathit{I}}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}-{\widehat{\lambda }}_k{\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}={\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}-{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}$，则

$$\begin{array}{c}{\widehat{\mathit{z}}}_{{\mathit{I}}_{\mathit{k}}}={\mathit{M}}_{\mathit{I}}{\widehat{\mathit{d}}}_{{\mathit{I}}_{\mathit{k}}}-{\widehat{\lambda }}_k{\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}+{\mathit{M}}_{\mathit{I},\mathit{J}}{\widehat{\mathit{d}}}_{{\mathit{J}}_{\mathit{k}}}+{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}=\\ {\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}-{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}+{\mathit{M}}_{\mathit{I},\mathit{J}}{\widehat{\mathit{d}}}_{{\mathit{J}}_{\mathit{k}}}+{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}=\\ {\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}-{\mathit{q}}_{{\mathit{I}}_{\mathit{k}}}-{\mathit{M}}_{\mathit{I},\mathit{J}}{\mathit{x}}_{{\mathit{J}}_{\mathit{k}}}+\mathit{q}{\mathit{I}}_{\mathit{k}}=\mathrm{0}\\ \end{array}$$

可知${\widehat{\mathit{z}}}_{\mathit{k}}=\left(\begin{array}{c}{\widehat{\mathit{z}}}_{{\mathit{J}}_{\mathit{k}}}\\ \mathrm{0}\end{array}\right)$，所以${\widehat{\mathit{z}}}_{\mathit{k}}^{\mathrm{T}}({\widehat{\mathit{d}}}_{\mathit{k}}+{\mathit{x}}_{\mathit{k}})=\mathrm{0}$。证毕。

定理4表明可以通过积极集算法找到问题（7）的解，即子问题（5）的解。积极集算法相当于是个组合问题，这保证了算法会在有限步结束，接下来讨论算法中一些具体问题，如在算法迭代步③中求解出${\beta }_k$使得$g(\mathit{d}({\beta }_k))=\underset{\mathrm{0}\le {\beta }_k\le {\beta }^{\mathrm{*}}}{\mathrm{m}\mathrm{i}\mathrm{n}}g(\mathit{d}(\beta ))$。

定理 5  　在积极集算法中，若$g(\mathit{d}({\beta }_k))=\underset{\mathrm{0}\le {\beta }_k\le {\beta }^{\mathrm{*}}}{\mathrm{m}\mathrm{i}\mathrm{n}}g(\mathit{d}(\beta ))$，则

$${\beta }_k=\left\{\begin{array}{c}{\beta }^{\mathrm{*}},{\beta }^{\mathrm{*}}\le {\beta }^{\mathrm{0}}\\ {\beta }^{\mathrm{0}},\mathrm{0}<{\beta }^{\mathrm{0}}<{\beta }^{\mathrm{*}}\\ \mathrm{0},{\beta }^{\mathrm{0}}\le \mathrm{0}\end{array}\right.$$

（14）

其中${\overline{\mathit{d}}}_{\mathit{k}}={\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}}$，${\beta }^{\mathrm{0}}=-\frac{{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}+{\mathit{q}}^{\mathrm{T}}{\overline{\mathit{d}}}_{\mathit{k}}}{{\overline{\mathit{d}}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}}$，${\beta }^{\mathrm{*}}=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\frac{({\mathit{x}}_{\mathit{k}}{)}_i+({\mathit{d}}_{\mathit{k}}{)}_i}{({\mathit{d}}_{\mathit{k}}{)}_i-({\widehat{\mathit{d}}}_{\mathit{k}}{)}_i}|i\in S_k\right\}$。

证明  　令${\overline{\mathit{d}}}_{\mathit{k}}={\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}}$，则$\mathit{d}(\beta )={\mathit{d}}_{\mathit{k}}+\beta ({\widehat{\mathit{d}}}_{\mathit{k}}-{\mathit{d}}_{\mathit{k}})={\mathit{d}}_{\mathit{k}}+\beta {\overline{\mathit{d}}}_{\mathit{k}}$。而

$$\begin{array}{c}g(\mathit{d}(\beta ))=g({\mathit{d}}_{\mathit{k}}+\beta {\overline{\mathit{d}}}_{\mathit{k}})=\\ \begin{array}{l}\frac{\mathrm{1}}{\mathrm{2}}({\mathit{d}}_{\mathit{k}}+\beta {\overline{\mathit{d}}}_{\mathit{k}}{)}^{\mathrm{T}}\mathit{M}({\mathit{d}}_{\mathit{k}}+\beta {\overline{\mathit{d}}}_{\mathit{k}})+\\ {\mathit{q}}^{\mathrm{T}}({\mathit{d}}_{\mathit{k}}+\beta {\overline{\mathit{d}}}_{\mathit{k}})=\frac{\mathrm{1}}{\mathrm{2}}({\overline{\mathit{d}}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}){\beta }^{\mathrm{2}}+\end{array}\\ ({\mathit{q}}^{\mathrm{T}}{\overline{\mathit{d}}}_{\mathit{k}}+{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}})\beta +\frac{\mathrm{1}}{\mathrm{2}}{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\mathit{d}}_{\mathit{k}}+{\mathit{q}}^{\mathrm{T}}{\mathit{d}}_{\mathit{k}}\end{array}$$

所以

$$\nabla g({\mathit{d}}_{\mathit{k}}+\beta {\overline{\mathit{d}}}_{\mathit{k}})=({\overline{\mathit{d}}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}})\beta +{\mathit{q}}^{\mathrm{T}}{\overline{\mathit{d}}}_{\mathit{k}}+{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}$$

又因为矩阵$\mathit{M}$是正定的，${\overline{\mathit{d}}}_{\mathit{k}}\ne \mathrm{0}$，则有${\overline{\mathit{d}}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}>\mathrm{0}$，所以$g(\mathit{d}(\beta ))$的最小值点为${\beta }^{\mathrm{0}}=$$-\frac{{\mathit{d}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}+{\mathit{q}}^{\mathrm{T}}{\overline{\mathit{d}}}_{\mathit{k}}}{{\overline{\mathit{d}}}_{\mathit{k}}^{\mathrm{T}}\mathit{M}{\overline{\mathit{d}}}_{\mathit{k}}}$。显然当${\beta }_k$满足式（14）时，使得$g(\mathit{d}({\beta }_k))=\underset{\mathrm{0}\le {\beta }_k\le {\beta }^{\mathrm{*}}}{\mathrm{m}\mathrm{i}\mathrm{n}}g(\mathit{d}(\beta ))$。证毕。

最后对问题（9）是否有解进行讨论。从引理1可以知道，当方程（10）有解时，显然问题（9）也有解。因为矩阵$\mathit{M}$是正定的，则它的主子矩阵${\mathit{M}}_{\mathrm{I}}$也是正定矩阵，即矩阵${\mathit{M}}_{\mathrm{I}}$可逆，那么有如下2种情况：①${\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}\ne \mathrm{0}$，则$\left|\begin{array}{cc}{\mathit{M}}_{\mathrm{I}}& {\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}\\ {\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}^{\mathrm{T}}& \mathrm{0}\end{array}\right|=\left|{\mathit{M}}_{\mathrm{I}}\right|\left|-{\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}^{\mathrm{T}}{\mathit{M}}_{\mathrm{I}}^{-\mathrm{1}}{\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}\right|\ne \mathrm{0}$。系数矩阵行列式不等于零，即矩阵可逆，方程有唯一解。②${\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}=\mathrm{0}$，则方程无解或有无穷多个解。

为了避免第2种情况出现，即保证${\mathit{b}}_{{\mathrm{I}}_{\mathit{k}}}\ne \mathrm{0}$，那么可对指标集$J_k$做如下处理。因为$\mathit{b}$为非零向量，则必可找到一个$i$使得${\mathit{b}}_{\mathit{i}}\ne \mathrm{0}$。若$i\in J_k$，则$J_k=J_k/i$；若$i\notin J_k$，则指标集$J_k$保持不变。这样保证了${\mathit{b}}_{{\mathit{I}}_{\mathit{k}}}\ne \mathrm{0}$，方程（10）有唯一解，所以在实际计算中可以保证问题（9）有解。

3 数值实验

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

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

积极集算法中求解子问题的初始向量${\mathit{d}}_{\mathrm{0}}$的公式如下：

$$({\mathit{d}}_{\mathrm{0}}{)}_i=\left\{\begin{array}{c}\frac{\frac{\mathrm{1}}{\mathrm{2}}{\mathit{x}}_{\mathrm{0}}^{\mathrm{T}}\mathit{B}{\mathit{x}}_{\mathrm{0}}+\mathrm{1}}{(\mathit{B}{\mathit{x}}_{\mathrm{0}}{)}_i}-({\mathit{x}}_{\mathrm{0}}{)}_i,j=i\\ -({\mathit{x}}_{\mathrm{0}}{)}_j,j\ne i\end{array}\right.$$

其中，$i$为$\mathit{B}{\mathit{x}}_{\mathrm{0}}$中第1个分量为正的指标。因为矩阵$\mathit{B}$为正定矩阵，这样的分量一定存在。这样选取保证了${\mathit{d}}_{\mathrm{0}}$在可行域内，也从侧面证明了子问题可行域非空。在第$k$次外迭代后，积极集算法中求解子问题的初始向量为${\mathit{d}}_{\mathit{k}}=(\mathrm{1}-{\alpha }_{k-\mathrm{1}}){\mathit{d}}_{\mathit{k}-\mathrm{1}}$，其中${\alpha }_{k-\mathrm{1}}$和${\mathit{d}}_{\mathit{k}-\mathrm{1}}$分别为上一步外迭代得到的步长和子问题的解。另外，积极集算法在迭代步③中（b）的情况2中$m$的选取也会影响迭代次数，在这里取$m=\left[\frac{l}{\mathrm{5}}\right]$，其中$l$表示${\widehat{z}}_k$小于零的分量的个数。

为了有效地分析算法SQP（J）和SQP（Q）性能，也用算法BAS^［16］求解了所有数值实验。在本节的所有数值实验中，初始矩阵$\mathit{A}$、$\mathit{B}$都为随机生成的对称正定矩阵。初始向量为${\mathit{x}}_{\mathrm{0}}={\mathit{e}}^{\mathit{s}}$，其中${\mathit{e}}^{\mathit{s}}$表示第$s$个分量为1、其余分量为零的向量，且$s=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\{r_i|i=\mathrm{1},\dots ,n\}$，$r_i=\mathrm{m}\mathrm{i}\mathrm{n}\{a_{ji}{b}_{ii}-a_{ii}{b}_{ji}|j=\mathrm{1},\dots ,n\}$。所有算法的停止准则为${‖{\mathit{d}}_{\mathit{k}}‖}_{\mathrm{2}}<{\mathrm{10}}^{-\mathrm{8}}$。另外，还讨论了算法的互补性，即$H={‖\mathrm{m}\mathrm{i}\mathrm{n}\{\mathit{A}{\mathit{x}}_{\mathit{k}}-{\lambda }_k\mathit{B}{\mathit{x}}_{\mathit{k}},{\mathit{x}}_{\mathit{k}}\}‖}_{\mathrm{2}}$。

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

$$\mathit{A}=\left[\begin{array}{cccc}\mathrm{4}& -\mathrm{2}& \mathrm{0}& \mathrm{0}\\ -\mathrm{2}& \mathrm{7}& \mathrm{6}& \mathrm{0}\\ \mathrm{0}& \mathrm{6}& \mathrm{9}& -\mathrm{1}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{8}\end{array}\right]$$，$$\mathit{B}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]$$

初始向量为${\mathit{x}}_{\mathrm{0}}={\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}$，通过算法SQP（J）和算法SQP（Q）得到

$${\mathit{x}}^{\mathit{*}}=\left[\begin{array}{c}\mathrm{1.264}\mathrm{9}\\ \mathrm{0.632}\mathrm{5}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right],$$$${\lambda }^{\mathrm{*}}=\mathrm{3}$$，$$\mathit{A}{\mathit{x}}^{\mathit{*}}-{\lambda }^{\mathrm{*}}\mathit{B}{\mathit{x}}^{\mathit{*}}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{3.794}\mathrm{7}\\ \mathrm{0}\end{array}\right]$$

这意味这${\lambda }^{\mathrm{*}}$和${\mathit{x}}^{\mathit{*}}$是矩阵对$\left(\mathit{A},\mathit{B}\right)$的Pareto特征值和对应的Pareto特征向量。

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

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

表1 随机生成的50阶问题的数值结果比较

Tab. 1 Comparison of numerical results of randomly generated 50th order problems

实验 次数	BAS			SQP（Q）			SQP（J）
	CPU运行时间/s	IT	H	CPU运行时间/s	IT	H	CPU运行时间/s	IT	H
1	0.000 947	65	4.445 8 ×10-6	0.022 173	8	2.428 9 ×10-6	0.006 272	8	1.947 9 ×10-6
2	0.001 478	126	6.103 4 ×10-6	0.033 262	8	3.979 4 ×10-6	0.006 002	8	2.73 76 ×10-7
3	0.000 890	71	4.753 9 ×10-6	0.023 004	8	1.101 5 ×10-5	0.006 316	8	8.552 7 ×10-7
4	0.003 668	109	6.513 7 ×10-6	0.021 788	8	1.1383 ×10-5	0.005 014	8	7.567 4 ×10-8
5	0.002 325	185	4.928 0 ×10-6	0.024 468	8	2.097 6 ×10-5	0.005 053	8	1.587 3 ×10-6

表2 随机生成的500阶问题的数值结果比较

Tab. 2 Comparison of numerical results of randomly generated 500th order problems

实验 次数	BAS			SQP（Q）			SQP（J）
	CPU运行时间/s	IT	H	CPU运行时间/s	IT	H	CPU运行时间/s	IT	H
1	0.647 598	775	4.505 1 ×10-5	0.328 313	10	6.479 5 ×10-4	0.243 969	10	1.719 2×10-9
2	0.765 895	903	4.902 7 ×10-5	0.337 919	10	8.700 5 ×10-4	0.250 600	10	1.308 9 ×10-9
3	0.825 336	969	5.971 3 ×10-5	0.347 302	10	4.131 4 ×10-4	0.227 094	10	1.936 7 ×10-9
4	0.092 072	100	4.926 5 ×10-6	0.366 384	10	9.921 4 ×10-5	0.298 834	10	1.745 1 ×10-9
5	0.265 748	298	2.343 4 ×10-5	0.358 240	10	2.519 9 ×10-4	0.232 598	10	1.223 7 ×10-9

表3 随机生成的1 000阶问题的数值结果比较

Tab. 3 Comparison of numerical results of randomly generated 1000th order problems

实验 次数	BAS			SQP（Q）			SQP（J）
	CPU运行时间/s	IT	H	CPU运行时间/s	IT	H	CPU运行时间/s	IT	H
1	6.687 324	1 360	8.917 8 ×10-5	1.444 580	11	0.002 1	1.764 132	10	2.272 9 ×10-7
2	0.503 070	103	4.808 7 ×10-6	1.460 458	10	0.001 2	1.650 242	10	2.788 6 ×10-7
3	5.928 741	1 197	8.351 4 ×10-5	1.426 617	10	0.002 5	1.347 867	10	1.959 7 ×10-7
4	8.955 474	1 812	1.141 1 ×10-4	1.442 109	10	0.001 1	1.861 311	10	2.946 4 ×10-7
5	6.787 376	1 353	9.138 6 ×10-5	1.468 115	10	0.002 1	1.632 466	10	1.945 3 ×10-7

4 结论

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

作者贡献声明

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

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

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