雷 淵，朱 琳，李 斌

（湖南大學(xué)數(shù)學(xué)學(xué)院，湖南長沙 410082）

互補(bǔ)特征值問題［1］（EiCP）是一類特殊的互補(bǔ)問題［2］，又稱為錐約束下的特征值問題［3-5］，它是經(jīng)典特征值問題的推廣。近年來，互補(bǔ)特征值問題受到了廣泛的關(guān)注，在結(jié)構(gòu)機(jī)械系統(tǒng)的動力學(xué)分析、信號處理、摩擦彈性系統(tǒng)和聲波系統(tǒng)的穩(wěn)定性分析以及流體動力學(xué)等方面都涉及到互補(bǔ)特征值問題的理論和對應(yīng)的特征值與特征向量的計算問題［6-8］。

目前關(guān)于互補(bǔ)特征值問題研究主要集中于2類常用的自對偶錐：非負(fù)錐和二階錐［9-10］。當(dāng)錐為非負(fù)錐Rn+時，該問題又稱為Pareto特征值問題，其數(shù)學(xué)描述為：給定n階實(shí)矩陣A和對稱正定矩陣B，求實(shí)數(shù)λ∈R和非零向量x∈Rn滿足

稱滿足式（1）的實(shí)數(shù)λ是矩陣對(A，B)的一個Pareto特征值，x是相應(yīng)的互補(bǔ)特征向量。當(dāng)式（1）中的矩陣A為實(shí)對稱矩陣時，該問題可稱為對稱Pareto特征值問題。由于Pareto特征值問題（1）等價于單純形上一個連續(xù)函數(shù)的變分不等式問題，所以該問題的解一定存在［11］。

與線性或非線性互補(bǔ)問題類似，利用非線性互補(bǔ)函數(shù)可以將Pareto特征值問題等價地轉(zhuǎn)化成一個非光滑方程組，然后用半光滑牛頓方法求解［12］，但是它在許多情況下會失敗。另外，求解約束規(guī)劃問題的內(nèi)點(diǎn)算法［13］、投影算法［14］以及DC（difference of convex functions）算法［15］都可用于計算大規(guī)模實(shí)對稱互補(bǔ)特征值，但都存在計算效率不高的問題。在文獻(xiàn)［16］中提出了一種具有譜選擇搜索方向的塊積極集算法（BAS），并通過大量的計算實(shí)驗(yàn)證明了其有效性。

序列二次規(guī)劃（SQP）算法是求解非線性規(guī)劃問題的一種高效算法［17-18］。序列二次規(guī)劃算法通過構(gòu)造一系列二次規(guī)劃子問題來求解非線性規(guī)劃問題，而積極集方法適合求解二次規(guī)劃問題。本文主要將積極集方法和序列二次規(guī)劃算法相結(jié)合，構(gòu)造求解對稱Pareto特征值問題的一類積極集方法。

1 求解對稱Pareto特征值的序列二次規(guī)劃方法

主要討論對稱Pareto特征值問題，即矩陣A、B都是對稱矩陣，其中B是對稱正定矩陣。當(dāng)矩陣A非正定，可以找到一個足夠大的正常數(shù)δ＞0，使得矩陣A+δB是正定的。若(λ，x)是矩陣對(A+δB，B)的Pareto特征值和特征向量，則由式（1）可以很容易地得到(λ?δ，x)是矩陣對(A，B)的一個Pareto特征值和特征向量。因此下面僅考慮A、B都是對稱正定矩陣的情況，用符號Sn+表示n階實(shí)對稱正定矩陣的集合。

在求解對稱Pareto特征值問題時，可以將原問題（1）轉(zhuǎn)化為如下約束優(yōu)化問題：

其拉格朗日函數(shù)為

顯然問題（2）的穩(wěn)定點(diǎn)滿足如下KKT（Karush-Kuhn-Tucker）條件：

式中：λ∈R、z∈Rn+為相應(yīng)等式約束和不等式約束的拉格朗日乘子。從式（4）可以知道(λ，x)是原問題（1）的一個解。因此可以利用序列二次規(guī)劃方法找到優(yōu)化問題（2）的一個穩(wěn)定點(diǎn)，繼而得到原問題（1）的一個Pareto特征值和其特征向量。

首先由當(dāng)前迭代點(diǎn)x k構(gòu)造二次規(guī)劃子問題：

其中矩陣M是個對稱正定矩陣。當(dāng)可行域非空時，問題（5）有解且有唯一解，記問題（5）的解為d k。用d k作為第k次迭代的搜索方向，再選取適當(dāng)步長αk，從而由迭代格式為x k+1=x k+αk d k生成序列{x k}。可以證明當(dāng)問題（5）的解d k=0時，x k就是原問題（1）的一個Pareto特征向量。

下面定義一個罰函數(shù)Pσ(x)：

式中：σ是罰參數(shù)，σ＞0。在經(jīng)典的序列二次規(guī)劃算法中，當(dāng)給定的罰參數(shù)σ對所有的k都滿足σ＞|λk|，其中λk是問題（5）等式約束的拉格朗日乘子，則d k是罰函數(shù)Pσ(x)的下降方向。又因?yàn)榱P函數(shù)Pσ(x)＞0有下界，故算法可收斂。在文獻(xiàn)［19］中有類似證明，這里不再闡述。在后面的數(shù)值實(shí)驗(yàn)中，發(fā)現(xiàn)這確實(shí)可以實(shí)現(xiàn)，其罰參數(shù)σ隨著每次迭代而更新。

簡單描述用SQP算法求解對稱正定矩陣Pareto特征值的基本框架如下。

（1）初始值。①給定對稱正定矩陣A，B∈Sn+；②選取對稱正定矩陣M∈Sn+和罰參數(shù)σ＞0；③選取誤差參數(shù)ε＞0和初始迭代點(diǎn)x0≥0，k=0。

（2）迭代。①求出子問題（5）唯一解d k；②如果，則算法停止，輸出x k。否則，轉(zhuǎn)至步驟③；③選取恰當(dāng)?shù)牟介Lαk∈(0，1]，使得Pσ(x k+αk d k)=轉(zhuǎn)至步驟①。

從上面算法可知，利用序列二次規(guī)劃問題求Pareto特征值問題需要找到適當(dāng)?shù)牟介Lαk∈(0，1]，使 得在 研 究過程中發(fā)現(xiàn)函數(shù)Pσ(x k+αk d k)在區(qū)間(0，1]上是分段連續(xù)函數(shù)，可以直接計算出需要的步長，這里就不贅述求解過程，直接給出結(jié)果（詳見文獻(xiàn)［19］）。

定理1設(shè)x k和d k由算法SQP生成，假設(shè)d k≠0且σ＞ρ，其中ρ是矩陣B?1A的譜半徑。如式（1）和式（2）選取的步長αk：

（1）當(dāng)ck＞0時

（2）當(dāng)ck≤0時

定理2設(shè){x k}是算法SQP生成的迭代序列，當(dāng)罰參數(shù)σ對所有的k都滿足σ＞max{|λk|，ρ}，其中λk是問題（5）等式約束的拉格朗日乘子，ρ是矩陣B?1A的譜半徑，則序列{x k}的任意聚點(diǎn)都是原問題（1）的一個Pareto特征向量。

從上面可以發(fā)現(xiàn)算法SQP的主要思想是用子問題（5）解d k作為搜索方向，生成迭代序列，最后收斂到原問題的解。如何快速精確求解子問題（5）是問題的關(guān)鍵，下面利用積極集方法求解子問題。

2 積極集算法

眾所周知，積極集方法被認(rèn)為是求解中小規(guī)模二次規(guī)劃問題的高效的算法之一。本節(jié)主要目的是構(gòu)造一類有效求解子問題（5）的解d k的積極集方法。已知當(dāng)前點(diǎn)x k，為了簡潔，將下面的量重新進(jìn)行標(biāo)記：，則可行域可簡記為D：={d|bTd=c，d+x k≥0}。進(jìn)而子問題（5）可以簡寫為

對于對稱正定矩陣M，當(dāng)可行域D非空時，易知問題（7）存在唯一解。那么d k是問題（7）的解當(dāng)且僅當(dāng)存在(λk，zk)滿足下面KKT條件：

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

下面定義2個指標(biāo)集：Jk=J(d k)={i|(d k)i+(x k)i=0，i∈N}，Ik=N/Jk，稱Jk為問題（7）在點(diǎn)d k處的積極集，其中N={1，2，…，n}。若d k是問題（7）的解，不難發(fā)現(xiàn)d k也是如下二次規(guī)劃問題的穩(wěn)定點(diǎn)：

引理1若問題（9）有解，則其解為

其中滿足下面方程：

證明不失一般性，對矩陣和向量進(jìn)行如下分塊：

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

則問題（9）轉(zhuǎn)化為

又因?yàn)閐 Jk=?x Jk，所以只需要求出的值即可。顯然可以通過解下面優(yōu)化問題得到要求的值。

問題（12）的解滿足KKT條件：

其中是其相應(yīng)的拉格朗日乘子，顯然解是下面方程的解：

證畢。

綜上可知，二次規(guī)劃問題（7）的解也是問題（9）的解，可通過求解問題（9）來找到問題（7）的解。由于問題（9）依賴于問題（7）的最優(yōu)解d k，所以解決問題（7）的關(guān)鍵步是確定哪些不等式約束是積極的，即設(shè)計一類確定指標(biāo)集Jk的策略。

下面建立積極集法求解子問題（7）的基本框架。

（1）初始值。①給定對稱正定矩陣M∈Sn+，向量x k∈Rn，q∈Rn，b∈Rn及標(biāo)量c∈R；②選取初始點(diǎn)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，轉(zhuǎn) 至 步驟①。

接下來討論積極集算法的一些基本性質(zhì)。

定理3設(shè)若序列{d k}由積極集算法生成，其中d0∈D，x k≥0，則有如下性質(zhì)：①對任意的k都有，d k∈D，即bTd k=c，d k+x k≥0；②對任意的k都有，g(d k+1)≤g(d k)，且等號成立當(dāng)且僅當(dāng)d k+1=d k。

證明（1）用歸納法證明。因?yàn)閐 0∈D，不妨假設(shè)d k∈D成立，所以只需證d k+1∈D即可。在積極集算法的情況（a）、（b）中d k+1=，d k+1是問題（9）的解，顯然有d k+1∈D。在情況（c）中分2種情況討論。

情況1：當(dāng)K k∩Sk=?時，因?yàn)閐 k+x k≥0，所以 對?i∈Sk，有d k+x k＞0。又 因 為Sk=所以對?i∈Sk，有?(x k)i＜(d k)i，進(jìn)而可知β*＞0。而

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

不等式約束分2種情況如下。①當(dāng)i?Sk時，易知。所以有

②當(dāng)i∈Sk時，易知。所以有

綜上可知d*+x k≥0，所以d*∈D。又因?yàn)閐 k+1=d(βk)，0≤βk≤β*，顯 然d k+1也 滿 足bTd k+1=c，d k+1+x k≥0，即d k+1∈D。

情況2：當(dāng)K k∩Sk≠?時，則d k+1是問題（13）的解，顯然有d k+1∈D。

（2）由上可知d k、d k+1都在可行域內(nèi)，由于d k和d k+1有多種取法，下面分情況討論。

①若d k+1=，那么d k+1是問題（9）的穩(wěn)定點(diǎn)，顯然d k也滿足問題（9）的約束條件，則易得g(d k+1)=即g(d k+1)≤g(d k)，等 號成立當(dāng)且僅當(dāng)d k+1=d k。②若d k+1=d(βk)=由 算 法 可 知g(d(βk))=顯然有g(shù)(d k+1)≤g(d k)，等號成立當(dāng)且僅當(dāng)d k+1=d k。③若d k+1=，即d k+1是問題（13）的解，顯然d k也滿足問題（13）的約束條件，同①可知g(d k+1)≤g(d k)，等號成立當(dāng)且僅當(dāng)d k+1=d k。

當(dāng)d k=且Jk=Jk?1/K k、或 者d k=d k+、或者這3種情況時，在這里就不一一討論了，與情況且=相同，可證g(d k+1)≤g(d k)，等號成立當(dāng)且僅當(dāng)d k+1=d k。證畢。

從上面定理3可知由積極集算法生成的序列{d k}在可行域內(nèi)，且保證目標(biāo)函數(shù)g(d)下降。

定理4若是算法輸出的解，則有如下性質(zhì)：，綜 上 可 知是問題（7）的解。

證明由積極集算法和定理3，性質(zhì)①②③④顯然易得，所以只需證明。因?yàn)?，?/p>

所以

定理4表明可以通過積極集算法找到問題（7）的解，即子問題（5）的解。積極集算法相當(dāng)于是個組合問題，這保證了算法會在有限步結(jié)束，接下來討論算法中一些具體問題，如在算法迭代步③中求解出βk使得

定理5在積極集算法中，若g(d(βk))=，則

證明令，則。而

所以

又因?yàn)榫仃嘙是正定的，則有所以g(d(β))的最小值點(diǎn)為。顯然當(dāng)βk滿足式（14）時，使得證畢。

最后對問題（9）是否有解進(jìn)行討論。從引理1可以知道，當(dāng)方程（10）有解時，顯然問題（9）也有解。因?yàn)榫仃嘙是正定的，則它的主子矩陣MI也是正定矩陣，即矩陣MI可逆，那么有如下2種情況：①bIk≠0，則。系數(shù)矩陣行列式不等于零，即矩陣可逆，方程有唯一解。②bIk=0，則方程無解或有無窮多個解。

為了避免第2種情況出現(xiàn)，即保證bIk≠0，那么可對指標(biāo)集Jk做如下處理。因?yàn)閎為非零向量，則必可找到一個i使得bi≠0。若i∈Jk，則Jk=Jk/i；若i?Jk，則指標(biāo)集Jk保持不變。這樣保證了b Ik≠0，方程（10）有唯一解，所以在實(shí)際計算中可以保證問題（9）有解。

3 數(shù)值實(shí)驗(yàn)

本節(jié)詳細(xì)描述數(shù)值算例的設(shè)計，嘗試通過數(shù)值實(shí)驗(yàn)驗(yàn)證積極集方法的有效性以及其在計算時間、迭代步數(shù)以及互補(bǔ)精度方面均優(yōu)于已有算法。所有數(shù)值實(shí)驗(yàn)在操作系統(tǒng)為64位windows 10和處理器為Intel（R）Core（TM）i7-8700 CPU@3.20GHz 3.19 GHz的電腦上進(jìn)行，使用的軟件是Matlab R2018b。

SQP算法最不可缺少的部分是子問題（5）中對稱正定矩陣M的選取。因?yàn)楸疚闹饕懻摼仃嘇是對稱正定的情況，在這里直接令M=A。為了簡潔，用積極集算法求解子問題的序列二次規(guī)劃問題，記為SQP（J）。又因?yàn)镸atlab軟件中的內(nèi)置函數(shù)“quadprog”能直接求解二次規(guī)劃，在所有實(shí)驗(yàn)中算法SQP（Q）表示直接用內(nèi)置函數(shù)“quadprog”求解子問題（5）。

積極集算法中求解子問題的初始向量d0的公式如下：

其中，i為Bx0中第1個分量為正的指標(biāo)。因?yàn)榫仃嘊為正定矩陣，這樣的分量一定存在。這樣選取保證了d0在可行域內(nèi)，也從側(cè)面證明了子問題可行域非空。在第k次外迭代后，積極集算法中求解子問題的初始向量為d k=(1?αk?1)d k?1，其中αk?1和d k?1分別為上一步外迭代得到的步長和子問題的解。另外，積極集算法在迭代步③中（b）的情況2中m的選取也會影響迭代次數(shù)，在這里取，其中l(wèi)表示小于零的分量的個數(shù)。

為了有效地分析算法SQP（J）和SQP（Q）性能，也用算法BAS［16］求解了所有數(shù)值實(shí)驗(yàn)。在本節(jié)的所有數(shù)值實(shí)驗(yàn)中，初始矩陣A、B都為隨機(jī)生成的對稱正定矩陣。初始向量為x 0=es，其中es表示第s個分量為1、其余分量為零的向量，且s=argmax{r i|i=1，…，n}，r i=min{ajibii?aiibji|j=1，…，n}。所 有算法的停止準(zhǔn)則為。另外，還討論了算法的互補(bǔ)性，即

例1為了測試算法的可行性，首先計算一個簡單的問題。在這里令

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

這意味這λ*和x*是矩陣對(A，B)的Pareto特征值和對應(yīng)的Pareto特征向量。

例2為了進(jìn)一步測試算法的有效性，測試一些隨機(jī)生成的問題。隨機(jī)生成5組階數(shù)分別為50、500和1 000的對稱正定矩陣對(A，B)，分別列出了算法BAS、算法SQP（Q）和算法SQP（J）計算矩陣對(A，B)的Pareto特征值所需的CPU運(yùn)行時間（s），迭代次數(shù)IT以及互補(bǔ)性H，見表1-3。其中，N表示矩陣的階數(shù)。

從表1、表2和表3可以看出在階數(shù)N=50時，算法BAS的運(yùn)行時間少于算法SQP（J），但互補(bǔ)性精度低于算法SQP（J）。當(dāng)矩陣階數(shù)增高到N=1000時，算法SQP（J）的運(yùn)行時間明顯低于算法BAS，而算法SQP（Q）的互補(bǔ)性表現(xiàn)糟糕。從整體運(yùn)算時間和互補(bǔ)性上來看，算法SQP（J）比算法BAS和算法SQP（Q）性能更穩(wěn)，效率更高。

表1 隨機(jī)生成的50階問題的數(shù)值結(jié)果比較Tab.1 Comparison of numerical results of randomly generated 50th order problems

表2 隨機(jī)生成的500階問題的數(shù)值結(jié)果比較Tab.2 Comparison of numerical results of randomly generated 500th order problems

表3 隨機(jī)生成的1 000階問題的數(shù)值結(jié)果比較Tab.3 Comparison of numerical results of randomly generated 1000th order problems

4 結(jié)論

構(gòu)造了求解實(shí)對稱互補(bǔ)特征值問題的積極集方法。實(shí)驗(yàn)數(shù)值結(jié)果表明，所提出的算法SQP（J）在迭代步數(shù)和計算時間以及互補(bǔ)性上均優(yōu)于算法BAS和算法SQP（Q）。算法SQP（J）能夠快速求解實(shí)對稱互補(bǔ)特征值問題，主要原因在于積極集算法更加適合求解子問題。如何選取子問題中的矩陣M是接下來需要研究的問題。

作者貢獻(xiàn)聲明：

雷 淵：提出研究思路，設(shè)計研究方案，修訂最終版本。

朱 琳：完成研究方案，進(jìn)行數(shù)值實(shí)驗(yàn)，分析數(shù)據(jù)，起草論文。

李 斌：負(fù)責(zé)數(shù)值實(shí)驗(yàn)的編程。