廣義嶺估計在矩陣病態(tài)問題中的應(yīng)用

王其騰，吳風(fēng)華

(華北理工大學(xué) 礦業(yè)工程學(xué)院，河北 唐山 063210)

測量數(shù)據(jù)的誤差處理是測量工作中最重要的環(huán)節(jié)之一，其對測量精度產(chǎn)生直接影響[1]。最小二乘(LS)估計在測量數(shù)據(jù)呈現(xiàn)線性回歸時，可以有效地對誤差進行消除[2]。但是若測量數(shù)據(jù)不是線性回歸，甚至呈現(xiàn)病態(tài)(即列向量之間存在嚴(yán)重的復(fù)共線型)時，LS估計處理的效果就很差[2-4]。故當(dāng)設(shè)計矩陣呈現(xiàn)病態(tài)時，需要設(shè)法消除其對最終測量精度的影響。有偏估計的提出就是為了消除矩陣病態(tài)影響，目前常用的有偏估計方法有嶺估計(Ridge Estimation)、廣義嶺估計(Generalized Ridge Estimate)、主成分估計(Principal Component Estimate)[5]。分別運用上述4種估計方法，利用MATLAB結(jié)合實際數(shù)據(jù)進行運算，驗證了廣義嶺估計在處理設(shè)計矩陣病態(tài)問題中有較另外2種具有偏估計更好的處理能力。

1設(shè)計矩陣病態(tài)及其判定

病態(tài)矩陣在測量數(shù)據(jù)處理中極其常見，因為在測量過程中，測量數(shù)據(jù)基本呈現(xiàn)非線性回歸，這樣就導(dǎo)致了設(shè)計矩陣各列向量之間可能存在嚴(yán)重的復(fù)共線性，稱之為病態(tài)矩陣[6]。當(dāng)矩陣呈現(xiàn)病態(tài)時，其所對應(yīng)的逆矩陣以及以該矩陣為系數(shù)的方程組的值都會因其微小的波動而產(chǎn)生較大的改變，對最終結(jié)果的運算帶來很大難度。

復(fù)共線性的診斷與度量是病態(tài)矩陣判斷的主要依據(jù)，目前比較常用的診斷方法有條件數(shù)法、特征分析法和方差膨脹因子法[7]。選擇條件數(shù)法作為診斷依據(jù)。

方N的條件數(shù)定義為：

(1)

式(1)中：λmax,λmin分別為N的最大及最小特征值。

根據(jù)上式，可很直觀地對N的特征值的離散程度進行判斷，從而進一步對復(fù)共線性是否存在以及矩陣病態(tài)程度進行判斷，最終對診斷模型病態(tài)性及病態(tài)程度進行判定。從以往大量的實際應(yīng)用生產(chǎn)中進行經(jīng)驗總結(jié)，可以大致地對條件數(shù)P進行如下判定：

(1)若0

(2)若100≤P≤1 000，則可以確定其存在較高程度的復(fù)共線性；

(3)若P>1 000，則其復(fù)共線性程度最強，病態(tài)最嚴(yán)重。

2原理與方案設(shè)計

2.1 原理基礎(chǔ)

現(xiàn)有線性平差模型：

(2)

式(2)中：設(shè)n為觀測值個數(shù)，N=BTB滿秩，同時t=rank(N),那么有X的最小二乘解[4]

(3)

具有唯一性，并且其為最優(yōu)的線性無偏估計解。但是在實際生產(chǎn)問題中，式(2)中N的行列式可能會很小，即det(N)≈0，數(shù)學(xué)中稱這種設(shè)計矩陣為病態(tài)矩陣，其所對應(yīng)的法方程為病態(tài)方程，其解非常不穩(wěn)定[8]。嶺估計、廣義嶺估計以及主成分估計的提出均是為了解決當(dāng)設(shè)計矩陣病態(tài)時，如何獲得更優(yōu)解的問題。

嶺估計是由最小二乘估計改進的一種有偏估計，其中參數(shù)X的嶺估計為[9]：

(4)

式(4)中：k>0為任意常數(shù)，稱其為嶺參數(shù)或偏參數(shù)。若k=0啊，則式(3)與式(4)相同。嶺估計相較于最小二乘估計最重要的改善是通過加入嶺參數(shù)K的概念來打破系數(shù)陣B的病態(tài)性，所以在進行嶺估計運算時，嶺參數(shù)K的確定就顯得尤其重要。許多統(tǒng)計學(xué)家提出了不同的確定嶺參數(shù)K的辦法，但是卻沒有一種辦法得到大家的公認，目前嶺跡法和雙h法是比較常用的方法[10]。該項研究選擇雙h法進行嶺參數(shù)K的確定，公式如下：

(5)

若將典則形式下的系數(shù)陣對角線上加上不同的參數(shù)k，將轉(zhuǎn)換為廣義嶺估計。故參數(shù)X的廣義嶺估計為：

(6)

上式中，G為對應(yīng)與N=BTB的特征根λ的標(biāo)準(zhǔn)正交化特征向量有：

GT(BTB)G=Λ=diag(λ1,…λt),K=diag(k,…kt)

(7)

式(7)中：k,…kt為t個嶺參數(shù)或稱廣義嶺參數(shù)。

主成分估計是將病態(tài)矩陣的一些趨近于0的特征根所對應(yīng)的成分從線性模型中剔除，對剩余的成分仍采用LS估計，最后根據(jù)關(guān)系X=GY確定X的估計。那么參數(shù)X的主成分估計為：

(8)

當(dāng)r=t時，式(8)與式(2)相同。

2.2 設(shè)計方案

圖1所示為數(shù)據(jù)處理設(shè)計方案。

圖1 數(shù)據(jù)處理設(shè)計方案

3算例計算

3.1 算例

在數(shù)據(jù)選擇過程中，因普通且未經(jīng)任何處理的數(shù)據(jù)并不能保證其具有病態(tài)效果，故選擇文獻[2]中所設(shè)計的方程系數(shù)陣，滿足其病態(tài)要求，同時因其為病態(tài)的方程系數(shù)矩陣，與實際生產(chǎn)過程中所產(chǎn)生的病態(tài)數(shù)據(jù)所對應(yīng)的系數(shù)陣基本相同，故其能代表普遍數(shù)據(jù)，滿足實驗要求。加入隨機誤差后的設(shè)計矩陣和觀測向量為[3,11]：

(9)

(9)

根據(jù)算例數(shù)據(jù)利用MATLAB編程，計算N=BTB的條件數(shù)為P=20 836.718 116 457 6，存在嚴(yán)重病態(tài)。遂利用上述3種有偏估計的程序分別進行運算。因選擇的數(shù)據(jù)為嚴(yán)重病態(tài)的數(shù)據(jù)，但為了展示效果，同時也運用最小二乘估計進行運算，運算結(jié)果只作為參考。

通過進行最小二乘估計以及上述3種有偏估計對算例數(shù)據(jù)進行運算，得到估計值如表1所示。

表1 不同估計方法估計結(jié)果

在實際工作中，無法通過參數(shù)估值來對估計方法消除矩陣病態(tài)的程度進行判定，所以對上述4種參數(shù)估值進行求均方誤差(MSE)和2-范數(shù)(F)運算，最后通過MSE及F對效果進行判定，下面對2-范數(shù)的定義進行簡單說明。

矩陣的2-范數(shù)，就是矩陣的轉(zhuǎn)置共軛矩陣與矩陣的積的最大特征根的平方根值，是指空間上2個向量矩陣的直線距離[12]。表2所示為不同估計方法MSE、2-范數(shù)結(jié)果。

分別將參數(shù)估計值、均方誤差及2-范數(shù)結(jié)果在圖中進行對比，如圖2、圖3所示。

圖2 4種估計方法估值對比

圖3 4種估計方法MSE、2-范數(shù)對比

3.2 算例分析

由于所選設(shè)計矩陣為嚴(yán)重病態(tài)，而最小二乘矩陣在處理矩陣病態(tài)問題中并不具備優(yōu)勢，故該研究中最小二乘估計僅作為參考，不參與最終的結(jié)果分析。

通過上述步驟可以發(fā)現(xiàn)，不同估計方法所得結(jié)果存在很大的差異，甚至在運用主成分估計進行運算時，x_3參數(shù)的符號也發(fā)生了變化，另外兩種有偏估計相對于最小二乘估計雖然并沒有發(fā)生符號的改變，但是其參數(shù)估計值也存在或大或小的差異。表明3種估計方法對矩陣病態(tài)的消除能力不盡相同。

在圖3中，通過對3種有偏估計以及最小二乘估計的均方誤差(MSE)和2-范數(shù)(F)進行分析，可以明顯地發(fā)現(xiàn)。在均方誤差(MSE)的對比中，嶺估計與最小二乘估計的消除病態(tài)能力大致相同，而主成分估計和廣義嶺估計有較大的消除能力，尤其是廣義嶺估計的消除效果更為明顯；在2-范數(shù)(F)的對比分析中，可以發(fā)現(xiàn)最小二乘估計、主成分估計、嶺估計與廣義嶺估計的消除病態(tài)能力逐步升高。

故通過選擇嚴(yán)重病態(tài)的算例數(shù)據(jù)，分別運用主成分估計(PCE)、嶺估計(RE)與廣義嶺估計(GRE)對算例中數(shù)據(jù)的處理，通過MSE及F-2直觀地分析出廣義嶺估計的消除矩陣病態(tài)的能力強于另外2種有偏估計方法。

4結(jié)論

(1)廣義嶺估計對典則形式下的系數(shù)陣對角線上加上了不同的廣義零參數(shù)k，這就導(dǎo)致廣義嶺估計的消除能力強于只添加嶺參數(shù)k的嶺估計，而主成分估計可能是因其去除了數(shù)據(jù)中非主要成分的影響，導(dǎo)致消除能力下降。

(2)當(dāng)數(shù)據(jù)呈現(xiàn)嚴(yán)重病態(tài)時，廣義嶺估計消除矩陣病態(tài)的能力是最強的，這就為處理矩陣病態(tài)數(shù)據(jù)時提供了更好的解決辦法。