王樂洋

1)東華理工大學(xué)測(cè)繪工程學(xué)院，南昌 330013

2)江西省數(shù)字國土重點(diǎn)實(shí)驗(yàn)室，撫州344000

1 引言

總體最小二乘(total least squares，TLS)自Golub 和Van Loan［1］于1980年首次從數(shù)值分析的觀點(diǎn)進(jìn)行分析并為之定名以來在算法和應(yīng)用方面得到了廣泛的研究［2－5］。在總體最小二乘解的性質(zhì)方面也有大量的研究，如總體最小二乘與最小二乘在解、殘差、數(shù)據(jù)擬合的改正數(shù)以及近似子空間的內(nèi)部聯(lián)系［2］，加權(quán)總體最小二乘問題的等價(jià)解集以及加權(quán)總體最小二乘解與加權(quán)最小二乘問題解之間的關(guān)系［6］，總體最小二乘解集、最小二乘解集以及與極小范數(shù)解之間的差異等等［7］。但是，關(guān)于總體最小二乘的擾動(dòng)分析研究不是很多，主要有總體最小二乘問題的統(tǒng)計(jì)特性和解擾動(dòng)的上限值的研究［2］;基于標(biāo)度總體最小二乘解存在且唯一的Golub－Van Loan 條件下的標(biāo)度總體最小二乘問題的擾動(dòng)分析［8］等。因?yàn)榭傮w最小二乘及標(biāo)度總體最小二乘問題并沒有最小二乘問題那樣簡單的范數(shù)表示的條件數(shù)，而分量條件數(shù)考慮了數(shù)據(jù)分量之間的關(guān)系，所以標(biāo)度總體最小二乘問題的擾動(dòng)分析是基于分量條件數(shù)的，而且由標(biāo)度總體最小二乘問題的條件數(shù)可以得到總體最小二乘問題的條件數(shù)［8］。從數(shù)值分析的角度進(jìn)行擾動(dòng)分析研究一般具有簡單、直接的優(yōu)點(diǎn)，但是，從測(cè)量平差和數(shù)值分析的角度進(jìn)行總體最小二乘、最小二乘及數(shù)據(jù)最小二乘解之間的擾動(dòng)分析及其關(guān)系的研究鮮有報(bào)道。標(biāo)度總體最小二乘擾動(dòng)分析的研究所給出的條件數(shù)是關(guān)于整個(gè)標(biāo)度總體最小二乘問題的，本文將從測(cè)量平差問題的法方程出發(fā)，從數(shù)值分析的角度，以系數(shù)矩陣條件數(shù)的定義為基礎(chǔ)并以奇異值分解為工具詳細(xì)推導(dǎo)總體最小二乘、最小二乘及數(shù)據(jù)最小二乘解的擾動(dòng)分析公式及三者之間的關(guān)系。

2 最小二乘的擾動(dòng)分析

2.1 普通最小二乘的擾動(dòng)分析

線性估計(jì)模型為

式中，A∈Rm×n(m ＞n)為列滿秩系數(shù)矩陣;X∈Rn×1為待估計(jì)參數(shù);b∈Rm×1為觀測(cè)值。

當(dāng)系數(shù)矩陣A 不含有誤差時(shí)，為LS 估計(jì)，法方程為

式中，N=ATA(非奇異)，L=ATb

式中，δL=AT(δb)。

當(dāng)NX=L≠0 時(shí)

所以

從而有

若矩陣N 為非奇異矩陣，則其條件數(shù)的定義為［9，10］

式中，‖·‖p表示矩陣的p 范數(shù)，p 取1、2 或∞。

式(8)說明:當(dāng)觀測(cè)值b 存在相對(duì)誤差(或擾動(dòng))時(shí)，將會(huì)引起LS 解的相對(duì)誤差，該相對(duì)誤差的上界是常數(shù)項(xiàng)‖δL‖/‖L‖的‖N‖‖N－1‖(即cond(N))倍，所以解的相對(duì)誤差的大小與條件數(shù)cond(N)的大小有關(guān);當(dāng)系數(shù)矩陣病態(tài)(cond(N)?1，即條件數(shù)相對(duì)較大)時(shí)，條件數(shù)cond(N)將對(duì)解的相對(duì)誤差起到放大作用。

2.2 數(shù)據(jù)最小二乘的擾動(dòng)分析

線性估計(jì)模型如式(1)所示，當(dāng)系數(shù)矩陣A 含有誤差而觀測(cè)值b 不含有誤差時(shí)的LS 估計(jì)稱為數(shù)據(jù)最小二乘(data least squares ，DLS)估計(jì)［11］，法方程如式(2)所示。設(shè)法方程系數(shù)矩陣N 的誤差(或擾動(dòng))為δN，式(1)的DLS 解為，系數(shù)矩陣N 的誤差(或擾動(dòng))δN 引起的解的擾動(dòng)為，模型精確解為X。

在上式兩端同乘以AT得

數(shù)據(jù)最小二乘的擾動(dòng)法方程為

即

所以

式中，δN=(δA)TδA，δA 為系數(shù)矩陣A 的誤差(或擾動(dòng))。

所以有

根據(jù)如下定理［10］:

如果‖B‖＜1，則I±B 為非奇異矩陣，且有估計(jì)

式中，‖·‖是矩陣的算子范數(shù)。

設(shè)‖N－1‖‖δN‖＜1，N+δN=N(I+N－1δN)為非奇異矩陣［9］，且有

由式(12)得

所以有

式(16)說明:如果條件數(shù)cond(N)(‖N‖‖N－1‖)越大，則系數(shù)矩陣N 的微小相對(duì)誤差‖δN‖/‖N‖將會(huì)引起解的相對(duì)誤差‖就越大，也就是對(duì)系數(shù)矩陣的相對(duì)誤差‖δN‖/‖N‖放大了‖N‖‖N－1‖倍。

2.3 觀測(cè)值與系數(shù)矩陣都存在誤差(或擾動(dòng))時(shí)的最小二乘擾動(dòng)分析

線性估計(jì)模型如式(1)所示，當(dāng)系數(shù)矩陣A 和觀測(cè)值b 同時(shí)含有誤差時(shí)，最小二乘(LS)估計(jì)的法方程為式(2)所示，設(shè)法方程系數(shù)矩陣N 的誤差(或擾動(dòng))為δN，觀測(cè)值b 的誤差(或擾動(dòng))為δb，δN 和δb 引起的解的擾動(dòng)為，模型精確解為X。

在上式兩端同乘以AT得

擾動(dòng)法方程為

即

所以

式中，δL=δATδb。

將式(17)展開結(jié)合式(2)得

對(duì)式(19)兩端取范數(shù)得

即

當(dāng)δN 足夠小時(shí)，有‖N－1‖‖δN‖＜1，則

由式(2)得

即

由式(22)和(24)得

即

由式(26)可以看出:當(dāng)系數(shù)矩陣A 和觀測(cè)值b同時(shí)含有誤差(或擾動(dòng))δA 和δb 時(shí)，最小二乘解LS的相對(duì)誤差是由兩部分引起的，即觀測(cè)值的相對(duì)誤差‖δL‖/‖L‖和法方程系數(shù)矩陣的相對(duì)誤差‖δN‖/‖N‖兩部分引起的;當(dāng)系數(shù)矩陣是病態(tài)矩陣時(shí)，條件數(shù)非常大，將會(huì)對(duì)觀測(cè)值的相對(duì)誤差‖δL‖/‖L‖和法方程系數(shù)矩陣的相對(duì)誤差‖δN‖/‖N‖起到放大作用，即引起LS 解的極大變化(或擾動(dòng))，也就是說‖δL‖/‖L‖和‖δN‖/‖N‖相同的條件下，條件數(shù)cond越大，LS 解的相對(duì)誤差也越大。

3 總體最小二乘的擾動(dòng)分析

線性估計(jì)模型如式(1)所示，當(dāng)系數(shù)矩陣A 和觀測(cè)值b 同時(shí)含有誤差時(shí)的估計(jì)為總體最小二乘(TLS)估計(jì)，法方程為［2］

式中，N=ATA(非奇異)，L=ATb，σn+1為增廣矩陣［A b］的最小奇異值，即有如下奇異值分解

式中，U=［U1U2］，U1=［u1，…un］，U2=［un+1…um］，ui∈Rm×1，UTU=Im;

設(shè)法方程系數(shù)矩陣N 的誤差(或擾動(dòng))為δN，觀測(cè)值b 的誤差(或擾動(dòng))為δb，式(1)的TLS 解為，δN 和δb 引起的解的擾動(dòng)為，模型精確解為X，擾動(dòng)法方程兩邊舍掉一次擾動(dòng)項(xiàng)ATδAX、ATδAδX、(δA)TAX、(δA)TAδX、2σn+1δσn+1X、2σn+1δσn+1δX、ATδb 和(δA)Tb，則

式中，δL=δATδb，δσn+1是增廣矩陣［δA δb］的最小奇異值。

將式(29)展開結(jié)合式(27)得

將式(31)兩端同乘以N－1得

將式(32)兩端取范數(shù)得

即

以矩陣的2-范數(shù)(譜范數(shù))進(jìn)行研究，則有［9］

式中，λATA=λN為矩陣ATA(即N)的最大特征值。

對(duì)系數(shù)矩陣A 進(jìn)行奇異值分解得［2，10］

式中，U'=［U'1U'2］，U'1=［u'1…u'n］，U'2=［u'n+1…u'm］，u'j∈Rm×1，U'TU'=Im;

V=［v'1… v'n］，v'i∈Rn×1，V'TV'=In;Σ'=diag(σ'1，…，σ'n)∈Rm×n;σ'1≥…≥σ'n≥0。

因?yàn)榫仃嘇 的非零奇異值是矩陣ATA(即N)的非零特征值的正平方根，所以結(jié)合式(36)得

根據(jù)奇異值交織定理得［2，12］

因此

假設(shè)擾動(dòng)足夠小，則有

所以

在式(43)兩邊同除以‖X‖得

由式(27)得

即

由式(44)和(46)得

由式(47)可以看出:當(dāng)系數(shù)矩陣A 和觀測(cè)值b同時(shí)含有誤差時(shí)的總體最小二乘(TLS)解的相對(duì)誤差是由三部分引起的，即觀測(cè)值的相對(duì)誤差‖δL‖/‖L‖、法方程系數(shù)矩陣的相對(duì)誤差‖δN‖/‖N‖以及系數(shù)矩陣和觀測(cè)值組成的增廣矩陣的擾動(dòng)當(dāng)系數(shù)矩陣病態(tài)時(shí)，條件數(shù)cond(N)(‖N‖‖N－1‖)非常大，將會(huì)對(duì)觀測(cè)值的相對(duì)誤差‖δL‖/‖L‖和法方程系數(shù)矩陣的相對(duì)誤差‖δN‖/‖N‖以及增廣矩陣的擾動(dòng)起到放大作用，即引起TLS 解的極大變化(或擾動(dòng))，也就是說在‖δL‖/‖L‖、‖δN‖/‖N‖和(δσn+1)2‖‖L‖相同的條件下，條件數(shù)cond(N)(‖N‖‖N－1‖)越大，TLS 解的相對(duì)誤差也越大。

同時(shí)，式(47)是前面各種最小二乘(LS)擾動(dòng)分析的概括(統(tǒng)一)形式，具體為:

1)當(dāng)系數(shù)矩陣不含誤差且沒有擾動(dòng)，僅觀測(cè)值含有誤差(或擾動(dòng))時(shí)，進(jìn)行LS 估計(jì)，‖δN‖=0，，則式(47)變?yōu)槭?8);

2)當(dāng)觀測(cè)值不含誤差且沒有擾動(dòng)，僅系數(shù)矩陣含有誤差(或擾動(dòng))時(shí)，進(jìn)行LS 估計(jì)(DLS 估計(jì))，‖δL‖，則式(47)變?yōu)槭?16);

3)當(dāng)觀測(cè)值和系數(shù)矩陣都含有誤差(或擾動(dòng))時(shí)，進(jìn)行LS 估計(jì)，，則式(47)變?yōu)槭?26)。

4 結(jié)束語

在實(shí)際參數(shù)估計(jì)問題中，一般來說數(shù)據(jù)采樣大小、模型化及測(cè)量等原因會(huì)引起系數(shù)矩陣的誤差，而總體最小二乘方法是可以同時(shí)顧及觀測(cè)值和系數(shù)矩陣誤差的有效的數(shù)據(jù)處理方法，因而在參數(shù)估計(jì)中得到廣泛的研究和應(yīng)用。本文從數(shù)值分析的角度出發(fā)，以系數(shù)矩陣條件數(shù)的定義為基礎(chǔ)并以奇異值分解為工具詳細(xì)推導(dǎo)了總體最小二乘、最小二乘及數(shù)據(jù)最小二乘解的擾動(dòng)分析公式，通過對(duì)比分析發(fā)現(xiàn)總體最小二乘擾動(dòng)分析公式是各種最小二乘擾動(dòng)分析公式的概括(統(tǒng)一)形式。

1 Golub G Hand and Van Loan C F.An analysis of the total least squares problem［J］.SIAM J Numer Anal.，1980，17:883－893.

2 Van Huffel S and Vandewalle J.The total least squares problem:Computational aspects and analysis［M］.SIAM，Philadelphia，1991.

3 Van Huffel S(Ed.).Recent advances in total least squares techniques and errors－in－variables modeling［M］.SIAM，Philadelphia，1997.

4 Van Huffel Sand Lemmerling P(Eds.).Total least squares and errors－in－variables modeling:Analysis，algorithms and applications［M］.Dordrecht:Kluwer Academic Publishers，2002.

5 王樂洋.總體最小二乘性質(zhì)研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，(5):48－52，57.(Wang Leyang.Research on properties of total least squares estimation［J］.Joural of Gesdesy and Geodynamics，2012;(5):48－52，57)

6 魏木生，陳果良.加權(quán)總體最小二乘問題的解集和性質(zhì)［J］.高校應(yīng)用數(shù)學(xué)學(xué)報(bào)A 輯，1994，9(3):304－311(Wei Musheng and Chen Guoliang.Solution sets and property for weighted total least squares problem［J］.Applied Mathematics－A Journal of Chinese Universities，1994，9(3):304－311.)

7 劉永輝，魏木生.TLS 和LS 問題的比較［J］.計(jì)算數(shù)學(xué)，2003，25(4):479－492(Liu Yonghui and Wei Musheng.On the comparision of the total least squares and the least squares problems［J］.Mathematica Numerica Sinica，2003，25(4):479－492.)

8 Zhou Liangmin，et al.Perturbation analysis and condition numbers of scaled total least squares problems，［J］Numer Algorithms:2009，51:381－399.

9 張賢達(dá).矩陣分析與應(yīng)用［M］.北京:清華大學(xué)出版社，2004(Zhang Xianda.Matrix analysis and applications［M］.Beijing:Tsinghua University Press，2004.)

10 易大義，沈云寶，李有法.計(jì)算方法(第二版)［M］.杭州:浙江大學(xué)出版社，2002(Yi Dayi，Shen Yunbao and Li Youfa.computational method(Second Edition)［M］.Hangzhou:Zhejiang University Press，2002.)

11 Paige C C and Strako? Z.Scaled total least squares fundamentals［J］.Numerische Mathematik，2002，91:117－146.

12 Thompson R C.Principal submatrices IX:Interlacing inequalities for singular values of submatrices［J］.Linear Algebraappl.，1972，5:1－12.