總體最小二乘的擾動分析*

王樂洋

1)東華理工大學測繪工程學院，南昌 330013

2)江西省數字國土重點實驗室，撫州344000

1 引言

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

2 最小二乘的擾動分析

2.1 普通最小二乘的擾動分析

線性估計模型為

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

當系數矩陣A 不含有誤差時，為LS 估計，法方程為

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

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

當NX=L≠0 時

所以

從而有

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

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

式(8)說明:當觀測值b 存在相對誤差(或擾動)時，將會引起LS 解的相對誤差，該相對誤差的上界是常數項‖δL‖/‖L‖的‖N‖‖N－1‖(即cond(N))倍，所以解的相對誤差的大小與條件數cond(N)的大小有關;當系數矩陣病態(cond(N)?1，即條件數相對較大)時，條件數cond(N)將對解的相對誤差起到放大作用。

2.2 數據最小二乘的擾動分析

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

在上式兩端同乘以AT得

數據最小二乘的擾動法方程為

即

所以

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

所以有

根據如下定理［10］:

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

式中，‖·‖是矩陣的算子范數。

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

由式(12)得

所以有

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

2.3 觀測值與系數矩陣都存在誤差(或擾動)時的最小二乘擾動分析

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

在上式兩端同乘以AT得

擾動法方程為

即

所以

式中，δL=δATδb。

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

對式(19)兩端取范數得

即

當δN 足夠小時，有‖N－1‖‖δN‖＜1，則

由式(2)得

即

由式(22)和(24)得

即

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

3 總體最小二乘的擾動分析

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

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

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

設法方程系數矩陣N 的誤差(或擾動)為δN，觀測值b 的誤差(或擾動)為δb，式(1)的TLS 解為，δN 和δb 引起的解的擾動為，模型精確解為X，擾動法方程兩邊舍掉一次擾動項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)展開結合式(27)得

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

將式(32)兩端取范數得

即

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

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

對系數矩陣A 進行奇異值分解得［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。

因為矩陣A 的非零奇異值是矩陣ATA(即N)的非零特征值的正平方根，所以結合式(36)得

根據奇異值交織定理得［2，12］

因此

假設擾動足夠小，則有

所以

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

由式(27)得

即

由式(44)和(46)得

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

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

1)當系數矩陣不含誤差且沒有擾動，僅觀測值含有誤差(或擾動)時，進行LS 估計，‖δN‖=0，，則式(47)變為式(8);

2)當觀測值不含誤差且沒有擾動，僅系數矩陣含有誤差(或擾動)時，進行LS 估計(DLS 估計)，‖δL‖，則式(47)變為式(16);

3)當觀測值和系數矩陣都含有誤差(或擾動)時，進行LS 估計，，則式(47)變為式(26)。

4 結束語

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

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 王樂洋.總體最小二乘性質研究［J］.大地測量與地球動力學，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 魏木生，陳果良.加權總體最小二乘問題的解集和性質［J］.高校應用數學學報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］.計算數學，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 張賢達.矩陣分析與應用［M］.北京:清華大學出版社，2004(Zhang Xianda.Matrix analysis and applications［M］.Beijing:Tsinghua University Press，2004.)

10 易大義，沈云寶，李有法.計算方法(第二版)［M］.杭州:浙江大學出版社，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.