帶狀結構的高斯圖模型嵌套懲罰估計

李凡群，張新生

（1.安徽財經大學 統計與應用數學學院，安徽 蚌埠 233030；2.復旦大學 管理學院，上海 200433）

0 引言

無向圖可以用G=(V，E)來描述，其中V包含p維隨機變量X=(X1，…，Xp)的p個頂點（node），邊集(edge)E={eij|1≤i＜j≤p} 描述了X=(X1，…，Xp)這p個頂點之間的兩兩條件相依性。當且僅當Xi和Xj在給定其他變量時是條件獨立的，則稱Xi和Xj之間不存在邊，即E中不存在元素eij。如果X=(X1，…，Xp)服從高斯分布N(μ，Σ)，則記 Ω=Σ-1為協方差矩陣的逆矩陣，稱為精確矩陣。在高斯分布場合的圖模型中，Xi和Xj之間不存在邊等價于給定其余變量的條件下，Xi和Xj是條件獨立的，這又等價于Ω中元素ωij=0。于是高斯圖模型的恢復和參數估計問題就等價于其協方差矩陣或其精確矩陣的估計以及精確矩陣中的零元素的識別問題。

Lauritzer[1]介紹了高斯圖模型的性質，并把圖模型選擇和估計的方法應用于高斯圖模型。Meinshausen和Bulmann[2]提出相鄰點選擇的方法。Levina等和Fredman等[3,4]把懲罰似然的方法用于精確矩陣的估計及圖模型的恢復。Fan和Peng[5]用自適應Lasso和SCAD懲罰選擇圖模型和參數估計。Zhang等(2014)[6]提出D-trace損失，在約束條件下，通過最小化D-trace損失，獲得了稀疏的精確矩陣的估計。

尤其是多維隨機變量有序的，比如縱向數據、時間序列、空間數據及光譜數據等。該類數據的特點是頂點間的條件相依性與頂點之間的位置有關，位置越遠，條件相依性越弱，稱該類模型為帶狀結構高斯圖模型。對這類數據，Smith和 Kohn[7]，Huang等[8]提出對正定矩陣做修正的Cholesky分解：Ω=(1-Φ)TD-2(1-Φ)，其中Φ是主對角線上的元素為1的下三角矩陣。D為對角形矩陣。Bickel等(2008)[9]通過對協方差矩陣的元素設置門限以確定帶寬，Levina等(2008)[3]對修正的Cholesky分解因子Φ提出嵌套的Lasso估計。Cai等(2011)[10]提出l1懲罰CLME方法，通過解線性規劃，獲得協方差矩陣的估計。本文利用精確矩陣的諸列(行)元素的系列線性回歸模型的解釋，對系列回歸模型的回歸系數施加嵌套的懲罰，得到精確矩陣的估計。與有關文獻相比，該方法更能直接獲得估計結果。隨機模擬結果表明，對精確矩陣施加嵌套的懲罰的估計提高了估計的精度。

1 精確矩陣的諸列元素的系列線性回歸模型的解釋

令X=(X1，…，Xp)T為一個p維隨機變量，且X～N(μ，Σ)。為了表達方便，本文用Σi，j表示移去Σ的第i行和第j列的的子矩陣，Σi，j表示Σ的第i行，其中第j個元素被移去，Xi表示移去變量X的第i個變量后的p-1維的隨機向量。如果對X，μ，Σ做如下分塊：

則由多元統計知，在給定X(2)=x(2)的條件分布為[11]：

特別地，取X(1)=Xi，X(2)=Xi

那么在給定其余變量Xi=xi的條件下，Xi的條件分布為：

如果把Xi看成響應變量，其余變量Xi看成協變量，考慮如下的線性回歸模型：

此時令β(j)=(β1，…，βj-1，βj+1…，βn)T，?i與Xi相互獨立，并且在均方意義下：

另一方面，若令Ω為多元正態分布X～N(μ，Σ)的精確矩陣，則由ΣΩ=I，有：

對每一個i(i=1，2，...，p)，得到精確矩陣Ω的諸對角線上的元素ωii以及其第i列上其余元素Ωi，i與回歸問題（1）中回歸系數β(j)及殘差方差之間有如下關系：

這表明精確矩陣Ω的第i列的元素Ωi，i與回歸系數成比例，于是Ωi，i中的零元素與β(j)中的零元素相對應，反之亦然。從而高斯圖模型的參數估計與模型選擇可以轉化為對回歸問題（1）的模型選擇及對殘差的方差的估計。基于估計

2 精確矩陣嵌套的l1懲罰估計

2.1 嵌套懲罰函數

根據前文的分析知，精確矩陣Ω中每行元素有相應的回歸問題中回歸系數的解釋。并且模型恢復與參數估計的問題轉化為對回歸系數β(j)和殘差方差的選擇和估計問題。

基于模型（1）的β(j)和的負的對數似然函數為：

(j)懲罰似然的估計為：

Levina等(2008)[3]基于改進的Cholesky分解，對cholsesky分解因子每行元素提出嵌套的懲罰，得到嵌套的分解因子估計。

下面本文基于帶狀模型對應的回歸模型，對其系數β(j))提出如下的嵌套懲罰函數：

該懲罰函數相當于加權的l1懲罰函數。在該懲罰下，一方面，當 k＜ j 時，若0。這意味著精確矩陣Ω的第j行元素有估計當l＞j時，若這意味著精確矩陣Ω的第j行元素有估計于是通過該懲罰，得到的圖模型估計保證了其帶狀結構特性。

考慮到 |βj-1|，|βj+1|的分母為1，故對 |βj-1|，|βj+1|兩個參數的懲罰在尺度上和其余參數不一致，所以對P1(|β(j)|)做如下兩種改進：

2.2 局部平方近似以及牛頓迭代算法

本文利用迭代算法，通過最小化：

得到β(j)和估計。

第一步：給定β(j)的當前估計值,計算：

對每一個j=1，2，…，p重復以上兩步，直到滿足一定的準則，從而得到諸β(j)和的估計，根據前文的分析，得到精確矩陣的估計。由于函數不是關于β(j)的凸函數，為了克服計算上的困難，本文對諸進行局部平方近似，得到了的顯式表達式。

一方面，對于懲罰函數的分子部分，令：

則上述三個懲罰函數經過局部平方近后得到如下形式：

其中C1為常數為對角形矩陣:

其中C2為常數，為對角形矩陣:

其中C3為常數，為對角形矩陣:

3 隨機模擬

下面通過幾個例子，分別計算估計的損失，比較不同方法的估計的精度。本文用K-L損失來度量估計的精度：把本文的方法得到的估計Modified nested Lasso簡記為MNL，并把結果同以下方法的估計進行比較，包括基于Cholesky分解的nested-Lasso(CNL)[3]，g-Lasso[4]，以及由經驗協方差矩陣得到的精確矩陣的估計記為Sample。考慮p=30，50，100，200的情況，對于每一個例子中的每一個p，從真實的模型中產生50組樣本數據，每組樣本數據各包含200個觀測值。分別利用BIC和CV準則，選取調節參數λ。

在進行迭代過程中，當參數的估計的絕對值|βk|＜10-5時，令 |βk|=0。由于模型1是稀疏的，所以本文報告了BIC準則和CV準則下的估計結果。對于每一個模型，由式（2）至式（4），本文分別計算了三種迭代估計的結果，列出其中估計的最小損失以及損失標準差。所有結果列于表1和表2所示。表1和表2表明，基于精確矩陣元素的回歸模型的改進的嵌套懲罰估計一致優于所列的其他方法。

表1 K-L損失下不同的估計方法對模型1的估計結果

表2 K-L 損失下不同的估計方法對模型2 的估計結果

（1）模型1(AR(3))：精確矩陣Ω=(ωij),其中ωij=0.5|i-j|，|i-j|≤3，i≠j,以及ωii=1,其他場合ωij=0。

（2）模型2（完全模型）：精確矩陣 Ω=(ωij),其中ωij=0.5|i-j|。

4 結束語

本文在具有帶狀結構的高斯圖模型場合下，利用高斯圖模型的精確矩陣的諸元素可以利用系列線性回歸模型進行解釋的性質，對諸線性回歸系數施加嵌套的l1懲罰，得到精確矩陣的改進的嵌套Lasso估計。通過局部平方近似及牛頓迭代算法，能很方便得到估計結果。與基于改進的Cholesky分解的嵌套Lasso估計相比，本文的方法能更直接得到精確矩陣的估計，減少了估計損失。數值模擬結果表明精確矩陣的嵌套Lasso估計提高了估計速度和估計精度。