求解多約束矩陣最小二乘問題的數值方法?

呂照美，劉新國

(中國海洋大學數學科學學院，山東 青島 266100)

1 研究內容

本文研究三約束矩陣最小二乘問題

(1)

X∈εpsd∩P∩B

其中

εpsd={X∈Rn×n；XT=X且λmin(X)≥ε},

B={X∈Rn×n；L≤X≤U},

〈A,B〉表示A與B的內積,(A,B)= tr(ATB)，tr(A)表示A的跡。記號A≤B表示Aij≤Bij，其中1≤i≤m，1≤j≤n。

對應的，

L=L1+L2+…+Lt,

U=U1+U2+…+Ut,

(2)

L，U有上述分裂，其中Li和Ui分別是L和U在Si上的投影，并記

(3)

上述問題出現于統計和數量經濟學等領域。統計學中的一個經典例子為：求一個到樣本協方差矩陣最近的對稱正定并且有一定結構的矩陣。這種有特殊結構的協方差矩陣[1]經常出現在物理、生物、心理學和社會科學模型。經濟學中的一個例子[2]是獲得效用函數問題，在這種情形下，擬合矩陣必須有界且對稱，同時它的最小特征值必須大于某一特定的正數。

由于Dykstra交替投影算法本身收斂較慢且內迭代中SPG方法非單調，使得整個算法收斂更慢。且文獻[12]中設置內外迭代次數上限均為5 000次，數值實驗表明，SPG達到迭代次數上限時存在沒有收斂或沒有達到精度要求的情況。文獻[12]中要求A,BT列滿秩并假定可行集非空以保證解的存在性與唯一性，但并未給出判斷可行集是否為空的方法以及在可行集上得到的點是否為最優解的條件。

在本文中,我們把Majorize-Minimize(以下簡稱MM)方法與Dykstra算法結合，提出求解式(1)的一種交替投影迭代方法。MM方法是一種單調收斂方法，Han[13]利用對偶性推導了Dykstra方法，并分析了其收斂性。

2 理論結果

2.1 問題的等價形式

令BP=B∩P，則

(4)

其中

(5)

這里，li和ui如式(3)定義，顯然式(4)等價于式(1)。

2.2 可解性

[tf(X)+(1-t)f(Y)]-f(tX+(1-t)Y)=

(6)

定理1令Ω=εpsd∩BP,設Ω≠φ，那么

(1)f(X)在Ω上有最小點。

證明 顯然，B為有界閉凸集，P為凸集，而由

λmin(tX+(1-t)Y)≤tλmin(X)+(1-t)λmin(Y)，

t∈[0,1]，XT=X，YT=Y，εpsd為閉凸集，如果Ω≠φ，則Ω為有界閉凸集，而f(X)為連續可微函數，因此f(X)在Ω上有最小點。

證明畢。

推論1如果Ω≠φ，而A和BT列滿秩，那么式(1)的解存在且唯一[12]。

下面分析Ω≠φ的條件。

記δi=λmin(Wi)，υi=λmax(Wi)則

{yTWiy:y∈Rn，‖y‖2=1}=[δi,υi]。

從而，如果

則

定理2如果

(7)

那么Ω≠φ。

2.3 最優性條件

定理3X*∈Ω為最小點的充要條件為

tr{(X-Y)TAT(AX,B-C)BT}≥0,?X∈Ω。

(8)

證明 如果X*為最小值點，則對?t∈[0,1]有

從而(8)式成立，否則取充分小，則得矛盾。

假設(8)式成立，則?X∈Ω，

證明畢。

注2式(8)不易驗證，下面給出一個充分條件。

推論2如果X*∈Ω，且

tr{(X-X*)AT(AX*B-C)BT}≥0，?X∈Bp，

(9)

則X*是原問題的最小點。

li≤αi≤ui(i=1,…,t)。

(10)

其中

bi=〈AGiB,C〉。

又

一方面，如果式(10)成立，則

(11)

另一方面，如果式(11)成立，則

有式(10)成立。

因此，推論2的等價形式為

(12)

則X*為原問題的最優解。

3 算法

MM算法是求解最優化問題的重要方法。為了求解問題(1)，設當前迭代點為X(k)，那么

2tr{BT(X-X(k))TAT(AX(k)B-C)}+

(13)

顯然，gk(X(k))=f(X(k))，因此有MM算法：

(14)

tr{BT(X-X(k))TAT(AX(k)B-C)}=

tr{(X-X(k))TAT(AX(k)B-C)BT}

且

2tr {(X-X(k))TAT(AX(k)B-C)BT} =

其中Ck與X無關，而

(15)

由此可知式(14)等價于

(16)

因此提出算法MM-Dykstra Method如下。

步驟2：計算

步驟3:計算

計算ck+1=‖X(k+1)-X(k)‖，如果ck+1

Boyle和Dykstra證明了Dykstra算法的收斂性[14]。

由式(14)，gk(X(k+1))≤gk(X(k))≡f(X(k))，且gk(X(k+1))≥f(X(k+1))，所以{f(X(k))}單調遞減。因此，在每次迭代僅需要求解兩個子問題。

其一是

(17)

其中Q是一個正交矩陣，

Λ(1)=diag(μ1,…,μs),

Λ(2)=diag(μs+1,..,μn),

μ1≥…≥μn≥ε>μs+1≥…≥μn,

則式(17)的解為

(18)

因此，

Pεpsd(X)=QΛ(i)QT。

(19)

另一個是

(20)

Z=Z1+…+Zt，

(21)

Zi與Gi有相同的結構，從而有

(22)

(23)

其解為

(24)

則

注3有多種取法。

由

其中，

aij=〈AGiB,AGjB〉，

4 初始點

由于MM算法是一種投影梯度法，一般而言收斂較慢。Lopez和Raydan[16]指出，Dykstra算法有時收斂也很慢，可以使用文獻[10]中的加速技巧，這里我們考慮選取特殊的來實現對算法的加速。

因此可以如下構造初始迭代點

(2)取

(Gi是Hadamard積)，?i滿足li≤?i≤ui且

|?i|=min{|αi|:li≤αi≤ui}。

5 數值實驗

在例1和例2中，L為零矩陣，U為所有元素均為3的矩陣，TOLin=1×10-7，TOLout=1×10-8，MM算法的最大迭代次數為5 000，P≡Rn×n，且在本文中

即在兩種算法比較的過程中所用的停止準則一致。

例1 對于文獻[12]中的實驗例子：ε=0，

B=I，

利用文獻[11]中算法4.1得到的結果與本文算法的結果對比見表1。

表1 MM_Dykstra算法與SPG_Dykstra算法的時間、誤差、殘差比較

其中在兩種算法中X*表達相同，均為

例2 取C=50×rand(m,p)，ε=0.1，此處以10階為例，計算結果見表2。

表2 MM_Dykstra算法與SPG_Dykstra算法10階矩陣的誤差與時間比較

從例1和例2兩個例子的對比中，我們能夠看到MM_Dykstra算法無論是在殘差還是在時間上都有了一定改善。

6 結語

本文提出了MM_Dykstra算法來求解多約束矩陣方程最小二乘問題，該方法將三約束問題轉換為雙約束問題，有效地減少了交替迭代次數。數值實驗結果表明，與SPG_Dykstra方法相比，本文算法具有更快的求解效果，結果精度也相對更高。