馮 巖， 宋 珊， 張子明， 徐常青

（蘇州科技大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，江蘇 蘇州215009）

協(xié)方差矩陣是多元統(tǒng)計(jì)學(xué)中的基本概念，它反映多變量之間的線性相關(guān)性及變量分布離散程度。 方差測(cè)量單個(gè)隨機(jī)變量的變化（比如一個(gè)群體中人的身高變化），而協(xié)方差則是衡量?jī)蓚€(gè)隨機(jī)變量的變化程度（如一個(gè)群體中的一個(gè)人的身高和體重的變化情況）。 Magnus J R 在1978 年提出了在干擾協(xié)方差矩陣中具有未知參數(shù)的GLS 模型的最大似然估計(jì)[1]，Meng X 等人在1991 使用EM 獲得漸近方差協(xié)方差矩陣[2]，Odell P L等人在1966 年提出了生成樣本協(xié)方差矩陣的數(shù)值程序[3]，Danaher P 等人在2012 年提出了用Lasso 方法對(duì)協(xié)方差矩陣的逆進(jìn)行估計(jì)[4]。 矩陣的理論發(fā)展也非常迅速，在19 世紀(jì)末，矩陣?yán)碚擉w系已基本形成。 到20 世紀(jì)，矩陣?yán)碚摰玫搅诉M(jìn)一步的發(fā)展。 目前，它已經(jīng)發(fā)展成為在物理學(xué)、控制論、機(jī)器人學(xué)、生物學(xué)、經(jīng)濟(jì)學(xué)等學(xué)科有大量應(yīng)用的數(shù)學(xué)分支。 而矩陣?yán)碚撛诮y(tǒng)計(jì)中主要用于估計(jì)線性模型的參數(shù)、對(duì)估計(jì)和模型進(jìn)行比較和利用矩陣的廣義逆研究隨機(jī)向量之間的關(guān)系等。 分塊矩陣是處理階數(shù)較高的矩陣時(shí)常采用的技巧，也是數(shù)學(xué)在多領(lǐng)域的研究工具。 對(duì)矩陣進(jìn)行適當(dāng)分塊，可使高階矩陣的運(yùn)算轉(zhuǎn)化為低階矩陣的運(yùn)算，同時(shí)也使原矩陣的結(jié)構(gòu)顯得簡(jiǎn)單而清晰，從而能夠大大簡(jiǎn)化運(yùn)算步驟，或給矩陣的理論推導(dǎo)帶來方便。 筆者主要運(yùn)用矩陣分塊理論對(duì)協(xié)方差矩陣進(jìn)行研究。

協(xié)方差矩陣反映隨機(jī)變量離散度及不同變量線性相關(guān)性的一類重要的二維統(tǒng)計(jì)參數(shù)，其對(duì)角元為相應(yīng)變量方差，而非對(duì)角元為相應(yīng)兩個(gè)變量間的協(xié)方差，它反映變量之間的線性相關(guān)性。 統(tǒng)計(jì)與概率論專家威廉·費(fèi)勒把協(xié)方差矩陣稱為隨機(jī)向量的方差。

協(xié)方差陣分為總體協(xié)方差陣和樣本協(xié)方差陣。 總體協(xié)方差矩陣（population covariance）一般記為Σ=（σij），它通常為未知參數(shù)陣，樣本協(xié)方差矩陣（sample covariance）也稱為樣本散度矩陣（dispersion matrix），一般記為SN（N 為樣本點(diǎn)個(gè)數(shù)）。 已知樣本情形下SN可計(jì)算，由此可對(duì)總體協(xié)方差矩陣進(jìn)行估計(jì)，即SN≈Σ。一般情況下，樣本點(diǎn)個(gè)數(shù)越大，近似效果越好。 但這種近似效果好壞不僅取決于樣本點(diǎn)個(gè)數(shù)，它還依賴于母體分布和取樣方式。 另一方面，很多情況下足夠大的N 不僅增加取樣成本，而且不現(xiàn)實(shí)。 如人類在觀察星云系亮度變化、海平面變化規(guī)律、地球生物物種變化與氣候環(huán)境的關(guān)系時(shí)，獲取的樣本觀察值個(gè)數(shù)有限，這時(shí)就需要通過其他方式（如母體概率分布假設(shè)和已知變量相關(guān)性等）來對(duì)總體協(xié)方差矩陣進(jìn)行估計(jì)。

假設(shè)有隨機(jī)變量x 和y，且對(duì)應(yīng)有N 個(gè)樣本點(diǎn)（x1，y1），（x2，y2），…，（xN，yN）。 那么隨機(jī)變量x 的樣本方差（散度）為

若x=y，則記var（x）=cov（x，x）=σx2=σ2（x）。 顯然σx2≥0，且有

現(xiàn)考慮n 維隨機(jī)向量x，y∈Rn。 則x 和y 對(duì)應(yīng)的樣本協(xié)方差為

利用協(xié)方差矩陣奇異值分解（SVD），可求出高維數(shù)據(jù)點(diǎn)分布的主成分及其在其不同主成分方向上的投影的方差。若隨機(jī)向量x 對(duì)應(yīng)的協(xié)方差矩陣為奇異，那么隨機(jī)變量x1，x2，…，xn具有一定的線性相關(guān)性。這時(shí)隨機(jī)變量出現(xiàn)冗余，可通過對(duì)應(yīng)協(xié)方差矩陣的分塊來研究其冗余性，從而簡(jiǎn)化對(duì)應(yīng)的回歸模型。

1 協(xié)方差矩陣及其性質(zhì)

設(shè)x1，x2，…，xn為n 個(gè)隨機(jī)變量，那么x=（x1，x2，…，xn）T為n 維隨機(jī)向量。 記A＝（aij）∈Rn×n為x 對(duì)應(yīng)的協(xié)方差矩陣，即

顯然AT=A（aij=aji，?（i，j）），即A 為實(shí)對(duì)稱矩陣。 不難證明A 為一個(gè)Gram 矩陣，即存在一組向量β1，β2，…，βn∈RK，使

其中

這里r=rank（A）為矩陣A 的秩，λi為A 的特征值，且A 的特征值均非負(fù)。 （3）式中矩陣U 的列向量（A 的特征向量）為相應(yīng)主成分（如U 的第i 個(gè)列向量ui為第i 個(gè)主成分），D 的對(duì)角元λi為相應(yīng)主成分變量的方差。 每個(gè)特征向量為定義主成分的原始變量保存一組回歸權(quán)重，主成分對(duì)應(yīng)的變量不相關(guān)。 下面考慮兩種情形：

（1）r=n。 此時(shí)在給定分布情形下，隨機(jī)向量x 的密度函數(shù)有定義；

（2）0＜r＜n。 此時(shí)令y＝UTx=（y1，y2，…，yn）T，那么有

定理1設(shè)n 維隨機(jī)向量x 的協(xié)方差矩陣A 有分解（3）式，其中U＝[u1，u2，…，un]為正交矩陣，Yk=UkTx，U＝[U1，U2]，Y1∈Rr，Y2∈Rn-r。 那么有

（1）隨機(jī)變量y1，y2，…，yr（yi=uiTx）不相關(guān)；

（2）對(duì)i=1，2，…，r，隨機(jī)變量yi的方差σi2=λi＞0，且σ12≥σ22≥…≥σr2；

（3）對(duì)i=r+1，r+2，…，n，有Pr（yi=uiTμ）=1，其中μ 為x 的期望E[x]。

證明注意到

故有

ei∈Rn為第i 個(gè)元為1 的坐標(biāo)向量。因此對(duì)i≠j，有cov（yi，yj）=0，從而隨機(jī)變量y1，y2，…，yr（yi=uiTx）不相關(guān)，故結(jié)論（1）成立。 對(duì)于結(jié)論（2），由（5）式，σi2=cov（yi，yi）=λi，又因?yàn)棣?≥λ2≥…≥λr＞0，所以σ12≥σ22≥…≥σr2，故結(jié)論（2）成立。 現(xiàn)證明結(jié)論（3）：對(duì)?i=r+1，r+2，…，n，有E[yi]=E[uiTx]=uiTE[x]=uiTμ，且由（5）式有cov（yi，yi）=0，所以，yi=uiTμ，即Pr（yi=uiTμ）=1。 證畢。

由定理1 可知，多維隨機(jī)變量的冗余性可以通過考察其對(duì)應(yīng)的協(xié)方差矩陣的奇異性來刻畫。 進(jìn)一步，通過其對(duì)應(yīng)的協(xié)方差矩陣的正交對(duì)角化找到這樣的正交變換，其變換得到的部分變量不相關(guān)。

2 正態(tài)分布中協(xié)方差矩陣的分塊

2.1 矩陣正態(tài)分布的定義

正態(tài)分布是統(tǒng)計(jì)學(xué)中最基礎(chǔ)和最重要的分布。 一元正態(tài)分布大約在200 多年前提出，而多元正態(tài)分布理論的提出和發(fā)展也經(jīng)歷了80 多年。 多元正態(tài)分布等價(jià)于一個(gè)多元隨機(jī)向量的分布。 1996 年，多變量分析統(tǒng)計(jì)專家Kollo T 等人[5]提出了隨機(jī)矩陣（包括隨機(jī)向量）的版本。 矩陣版本是向量版本的“雙線性”擴(kuò)展，并且多變量結(jié)構(gòu)從協(xié)方差結(jié)構(gòu)獲得，該協(xié)方差結(jié)構(gòu)將被呈現(xiàn)為兩個(gè)方差矩陣的Kronecker 積。 然而，另一方面，通過選擇特定的協(xié)方差結(jié)構(gòu)，總是可以從多元正態(tài)分布中獲得矩陣正態(tài)分布。

多元正態(tài)分布至少有三種定義和刻畫方法，即利用密度函數(shù)、特征函數(shù)和通過分布特征（如矩函數(shù)等）。文中的方法將依賴于強(qiáng)調(diào)正態(tài)分布和線性（多線性）變換之間聯(lián)系的特征，也可以使用其他特征。

眾所周知，一元標(biāo)準(zhǔn)正態(tài)分布的密度函數(shù)定義為

并表示為U～N（0，1）。 由此得出E[U]=0 且D[U]=1。 一般的，若隨機(jī)變量X 服從均值μ、方差σ2的正態(tài)分布，那么變量U=（X-μ）/σ 服從標(biāo)準(zhǔn)正態(tài)分布，因此，X 與隨機(jī)變量

具相同分布，故其密度函數(shù)為

記為X～N（μ，σ2），從（7）式可以得出，kX 與kμ+kσU 具有相同的分布。 因此，kX～N（kμ，（kσ）2）。 現(xiàn)在令u=（U1，…，Up）T是一個(gè)由p 個(gè)獨(dú)立的同分布（i.i.d.）N（0，1）元素組成的向量。由于獨(dú)立性，從（6）式可以得到u的密度函數(shù)為

記為u～Np（0，I）。

現(xiàn)考慮X 為p 維向量，并記均值E[X]=μ，協(xié)方差矩陣為D[X]=Σ，其中Σ 為p 階半正定矩陣（記Σ≥0）。若Σ＞0（正定），則存在滿秩陣τ∈Rp×p，Σ=ττT。 不難證明X～Np（μ，Σ）當(dāng)且僅當(dāng)X 與下式定義的隨機(jī)變量具有相同分布

其中u～Np（0，I）。 若Σ＞0，則其密度函數(shù)為

現(xiàn)在引入了矩陣正態(tài)分布。

定義1設(shè)矩陣X=（Xij）∈Rm×n為隨機(jī)矩陣，記服從標(biāo)準(zhǔn)正態(tài)分布，記X～Nm，n（0，Im，In），若：（1）X·1，X·2，…，X·ni.i.d.Nm（0，Im）；（2）X1·，X2·，…，Xm·i.i.d.Nn（0，In）。即X 的行向量和列向量均服從標(biāo)準(zhǔn)正態(tài)分布。

定義2令Σ=ττ，Ψ=γγT，其中，τ：p×r，γ：n×s，矩陣X：p×n 被稱為關(guān)于參數(shù)（μ，Σ，Ψ）的矩陣正態(tài)分布，如果它和下式定義的隨機(jī)矩陣有相同的分布

其中μ：p×n 是非隨機(jī)常數(shù)矩陣，U：r×s 是s 由個(gè)i.i.d.于Nr（0，I）的向量Ui（i=1，2，…，s）組成。 如果X：p×n 是矩陣正態(tài)分布，則將其表示為X～Np，n（μ，Σ，Ψ）。 如果Σ，Ψ 是正定的，則（12）式中的τ 和γ 均是非奇異矩陣。以下討論排除情況Σ＝0，Ψ＝0。

推論1若X=（Xij）∈Rm×n，且X～Np，n（μ，Σ，Ψ），則有

（1）X·1，X·2，…，X·ni.i.d.Nm（μ·j，ΨjjΣ）；（2）X1·，X2·，…，Xm·i.i.d.Nn（μi·，ΣiiΨ）。

定理2令X～Np，n（μ，Σ，Ψ），對(duì)任意的A：q×p 和B：m×n，有

定理3X～Nm，n（0，Im，In）當(dāng)且僅當(dāng)vecX～Nmn（0，Imn）。

證明令，其中，X·j=（X1j，X2j，…，Xmj）T，j=1，2，…，n

所以

X～Nm，n（0，Im，In）當(dāng)且僅當(dāng)vecX～Nmn（0，Imn）。

推論2X～Np，n（μ，Σ，Ψ）當(dāng)且僅當(dāng)vecX～Npn（vecμ，Ψ?Σ）。

性質(zhì)1若X～Np，n（μ，Σ，Ψ），則vecX 密度函數(shù)為

證明令x=vecX。 因若X～Np，n（μ，Σ，Ψ），則vecX～Npn（vecμ，Ψ?Σ），因此，由式（11）知

又因?yàn)関ec（AXB）=（BT?A）vecX，且tr（AB）=tr（BA），所以

又因?yàn)閨Ψ?Σ|=|Ψ|m|Σ|n，所以

推論3若X～Nm，n（0，Im，In），則vecX 的密度函數(shù)為

2.2 矩陣正態(tài)分布的一些性質(zhì)

當(dāng)協(xié)方差矩陣為奇異矩陣時(shí)，原始隨機(jī)變量之間一定存在相關(guān)性，這時(shí)隨機(jī)變量存在冗余，需要考慮協(xié)方差矩陣的分塊，希望剔除相關(guān)變量，并找出不相關(guān)隨機(jī)變量。

考慮若X=（Xij）∈Rm×n的分塊，使用以下記號(hào)

推論4令X～Np，n（μ，Σ，Ψ），其中X，μ，Σ 和Ψ 如（15-19）式所定義，那么

（1）X11～Nr，s（μ11，Σ11，Ψ11）；（2）X12～Nr，n-s（μ12，Σ11，Ψ22）；（3）X21～Np-r，s（μ21，Σ22，Ψ11）；（4）X22～Np-r，n-s（μ22，Σ22，Ψ22）；（5）X·1～Np，s（μ·1，Σ，Ψ11）；（6）X·2～Np，n-s（μ·2，Σ，Ψ22）；（7）X1·～Nr，n（μ1·，Σ11，Ψ）；（8）X2·～Np-r，n（μ2·，Σ22，Ψ）。

證明（1）令A(yù)=（Ir，0），則

由定理2 知AμB＝μ11，AΣAT=Σ11，BTΨB=Ψ11，因此，

同理可證（6）、（7）、（8）。

推論4 說明：一個(gè)服從矩陣正態(tài)分布的隨機(jī)矩陣，其對(duì)應(yīng)的每個(gè)列（行）隨機(jī)向量、子矩陣也服從矩陣正態(tài)分布。 下面將進(jìn)一步證明：服從矩陣正態(tài)分布的隨機(jī)矩陣，其對(duì)應(yīng)的分塊條件分布也服從矩陣正態(tài)分布。

3 條件正態(tài)分布協(xié)方差矩陣的分塊

這里仍采用（15-19）式中的記號(hào)，那么有

且

下面來研究隨機(jī)向量x1與x2的相關(guān)性。 有結(jié)論：

定理4[6]若Σ22非奇異，則x1|x2～Nn1（μ1.2，Σ11.2）， 其中

Σ11.2稱為Σ 關(guān)于Σ22的Schur 補(bǔ)。

證明令

由性質(zhì)知Ax～Nn（Aμ，AΣAT），將代入得

由上式可得x1-Bx2╧x2。 所以有

其中

假設(shè)Σ11是非奇異的并注意到任何C，有

假設(shè)Σ22是奇異的，則存在正交矩陣和非奇異矩陣D，使得

那么

又x2～Nn（μ2，Σ22）?x2-μ2～Nn（0，Σ22），H2T（x2-μ2）～Nn（0，H2TΣ22H2）。

因?yàn)镠 為正交矩陣，所以

因此，H1TH1=In1，H2TH2=In2，H1TH2=H2TH1=0，所以

因此，可以得到

由定理4 可知：服從矩陣正態(tài)分布的隨機(jī)矩陣，其對(duì)應(yīng)的分塊條件分布也服從矩陣正態(tài)分布。

4 奇異協(xié)方差矩陣對(duì)應(yīng)的線性回歸模型

給定隨機(jī)向量x=（x1，x2，…，xn）T∈Rn對(duì)應(yīng)的協(xié)方差陣Σ。若Σ 可逆，則在給定x 服從某分布的條件下，其對(duì)應(yīng)的分布密度函數(shù)可以經(jīng)其分布參數(shù)（μ，Σ）表示。 但在Σ 奇異的條件下，其對(duì)應(yīng)的密度函數(shù)無法表示，該節(jié)討論Σ 奇異時(shí)對(duì)應(yīng)的密度函數(shù)[7-12]。

設(shè)x～N（μ，Σ），則存在正交矩陣U，使得Σ=UDUT，其中D=diag（λ1，…，λn）。

若λ1≥λ2≥…≥λn＞0，即Σ 正定（Σ 可逆），則有

更一般的，有

即Σ 奇異，rank（Σ）=r

（1）廣義逆來代替逆。 用Σ+表示Σ 的廣義逆，則

其中U1=（u1，…，ur）∈Rn×r，因此

（2）因?yàn)棣?UDUT，所以其中，D1=diag（λ1，…，λr）。

令y=UTx=（y1，y2，…，yn）T，則

所以?i=1，…，r，有σ2（yi）=λi；?i=r+1，…，n，有σ2（yi）=0。

因此，y1，y2，…，yr為兩兩互相獨(dú)立的隨機(jī)變量，yr+1，…，yn為常量。 令Y1=（y1，y2，…，yr）

5 協(xié)方差矩陣在多元線性回歸中的應(yīng)用

該節(jié)將給出幾個(gè)實(shí)例，以說明協(xié)方差矩陣在多元線性回歸模型刻畫及其參數(shù)估計(jì)中的應(yīng)用。

例1假設(shè)隨機(jī)向量x∈Rm的均值為μ，方差為Ω，且Ω 為正定的，這時(shí)可以定義一個(gè)關(guān)于x 的線性變換，使得變換后的隨機(jī)向量是標(biāo)準(zhǔn)的，即它的均值為0，協(xié)方差矩陣為Im。 令Ω1/2是任意的滿足Ω=Ω1/2（Ω1/2）T的矩陣，并令z=Ω-1/2（x-μ），其中Ω-1/2=（Ω1/2）-1，可以發(fā)現(xiàn)

和

例2對(duì)于一個(gè)多元線性回歸模型

其中y=（y1，y2，…，yN）T∈RN，X∈RN×r， β=（β1，β2，…，βN）T∈RN， ε=（ε1，ε2，…，εN）T∈RN，并且ε～NN（0，σ2In），用來表示對(duì)參數(shù)β 的一個(gè)估計(jì)，則則

對(duì)f（β）關(guān)于β 求導(dǎo)得f′（β）=2（XTX）β-2XTy，令f′（β）=0，得到（XTX）β=XTy，因此，可以推出=（XTX）-1XTy。

當(dāng)ε～NN（0，σ2C）時(shí)，其中C=diag（c12，c22，…，cN2）為一個(gè)對(duì)角矩陣，ci是常數(shù)，則

現(xiàn)在，定義一個(gè)新的線性模型

等效地y*=X*β+ε*，其中

那么ε*的協(xié)方差矩陣為

參數(shù)β 的估計(jì)值為

例3在例2 中得到了在多元回歸規(guī)模型y=Xβ+ε 中對(duì)β 的一個(gè)加權(quán)最小二乘估計(jì)，現(xiàn)在考慮廣義最小二乘估計(jì)。

當(dāng)ε～NN（0，σ2C）時(shí)，其中C∈RN×N是一個(gè)已知的正定矩陣，現(xiàn)在想要得到一個(gè)使得隨機(jī)誤差的協(xié)方差為σ2In的變換，可令T∈RN×N為滿足TTT=C 的矩陣，定義一個(gè)新的多元線性模型

等效地y*=X*β+ε*，其中y*=T-1y，X*=T-1X，ε*=T-1ε，那么ε*的協(xié)方差矩陣為

因此，參數(shù)β 的估計(jì)值為

6 結(jié)語

筆者主要研究了矩陣分塊方法在協(xié)方差矩陣中的應(yīng)用。 首先，論文介紹了協(xié)方差矩陣的基本性質(zhì)及正態(tài)分布隨機(jī)變量的協(xié)方差矩陣的分塊特征；然后，通過矩陣分塊方法，結(jié)合矩陣分解和矩陣廣義逆，介紹了奇異協(xié)方差矩陣對(duì)應(yīng)的線性回歸模型和矩陣分塊方法在多元線性回歸中的應(yīng)用。