基于SVD求解病態線性方程組的正則因子分步選取方法

汪強強， 劉海飛， 柳建新， 李 星， 施昕祎

(1.中南大學 地球科學與信息物理學院，長沙 410083；2.有色資源與地質災害探查湖南省重點實驗室，長沙 410083)

0 引言

地球物理反演問題大多數是非線性多元函數的極值問題，由于涉及高階偏導數矩陣的計算和存儲，采用直接法對其進行求解尤為困難，必須通過線性化手段將非線性問題轉換為解方程組的問題[1]。將非線性問題進行泰勒級數展開并略去二次以上的項可以得到線性方程組，但得到的線性方程組往往是病態的[2]。利用奇異值分解法先將病態線性方程組進行分解，然后采用正則化技術削弱小奇異值的影響，可以使求得的解更好地逼近真實解，能夠達到有效求解病態方程組的目的[3]。目前削弱小奇異值影響的正則化技術主要有:Wiggins零化小奇異值法和Tikhonov濾波正則化方法[4]，基于Tikhonov濾波正則化方法求取病態線性方程組的正則解，尋找最佳正則化因子成為奇異值分解法成功求解病態線性方程組的關鍵。Hansen[5]提出了L-曲線的概念，并證明了點(ln‖b-Axλ‖2,ln‖xλ‖2)所構成的L-曲線上曲率最大的點對應最佳正則化因子的位置；O’ Leary等[6]利用三次樣條平滑的L-曲線計算曲率從而確定最佳角點；J.Longina Castellanos等[7]采用三角形法確定L-曲線的最佳角點；Hansen等[8]提出了一種自適應修剪算法，確定L-曲線最佳角點的效果較好。筆者在前人工作的基礎上，給出了一種病態線性方程組正則化因子的分步選取方法，即先根據L-曲線的局部特征找到可能的候選角點，再從候選角點中找到最佳角點，進而確定最佳正則化因子，再將其代入濾波正則化方程即得病態線性方程組的全局最優解。最后，通過希爾伯特方程組對算法進行測試，驗證了本文算法的有效性。

1 奇異值分解法

地球反演中常因線性反演方程的嚴重病態特征導致計算結果存在較大誤差，奇異值分解法是求解病態線性方程組主要方法之一。對于線性方程組

Ax=b

(1)

(2)

當λ大于某數之后，若A+λI正定，則可以得到原問題的直接解而無需奇異值分解。但是λ的取值難以確定，而奇異值分解無需判定矩陣A是否對稱正定，對線性方程組的系數矩陣進行奇異值分解后，可以通過分解的結果得到不同λ的正則解。A的奇異值分解為

(3)

其中：左矩陣U=(u1,u2,…,un)，且U∈Rm×n，UTU=In,In為單位矩陣；右矩陣V=(v1,v2,…,vn)，且V∈Rn×n，VTV=In；對角陣S=diag(σ1,σ2,…,σn)，且奇異值σi以σ1≥σ2≥…≥σn≥0依次遞減。則方程(1)的解x為

(4)

其中：A+為系數矩陣A的廣義逆；S+為對角陣S的廣義逆。從式(4)可以看出，隨著i的增加奇異值σi逐漸趨于零時或者說矩陣A的條件數σ1/σn很大時，直接求解Ax=b很難求得精確解，采用正則化削弱奇異值的影響，可以使所得解更接近精確解。經正則化后，線性方程組的解為[9]式(5)。

(5)

選擇合適的正則化因子λ，可以壓制較小奇異值的影響。

圖1 L-曲線示意圖Fig.1 Schematic diagram of L-curve

2 基于L-曲線尋找最佳正則化因子

2.1 L-曲線

所謂L-曲線是指點(ln‖b-Axλ‖2,ln‖xλ‖2)在直角坐標系中構成的曲線圖，其中‖xλ‖2為正則化解的范數，‖b-Axλ‖2為正則化解的殘差范數。根據Tikhonov正則化方法，當正則化因子λ增大時，正則化解的范數‖xλ‖2減小，相應的殘差范數‖b-Axλ‖2增大，反之正則化解的范數‖xλ‖2增大，殘差范數‖b-Axλ‖2減小，如圖1所示。當λ取到較為合理的值時，正則化解的范數‖xλ‖2和殘差范數‖b-Axλ‖2能夠達到較好的平衡，這個正則化因子就是最佳正則化因子，它對應L-曲線的最佳角點。

基于L-曲線選取最佳角點的方法主要有最大曲率法、三角形法等，O’Leary等[6]提出，如果平滑連續的L-曲線由(ρ,η)表示，ρ和η分別表示解的范數和對應的殘差范數，并且ρ和η為二階可導的函數，那么可以通過式(6)直接計算曲率k(λ)；對于由許多離散點組成的L-曲線，先擬合一條關于離散點的三次樣條曲線，再通過曲率計算公式計算曲率，從而確定最大曲率點，最大曲率點對應L-曲線的最佳角點。

(6)

(7)

最佳角點對應定向面積為負且角度小于7π/8的頂點。

盡管最大曲率法和三角形法都可以基于L-曲線選定最佳角點，但它們各自都存在不足之處。對于由離散點構成的L-曲線，最大曲率法通過擬合一條關于離散點的三次樣條曲線來計算曲率，但這條樣條曲線對離散點的分布很敏感，并且有時不能完全考慮所有的離散點，導致計算得到的最大曲率點并非對應最佳正則化因子。三角形方法在選定最佳角點時，位于L-曲線垂直分支上的離散點都參與了計算，導致計算量大，復雜性高。針對以上問題，我們提出了一種正則化因子的分步選取方法，考慮L-曲線的斜率突變區可能存在局部角點這一問題，基于一系列子L-曲線選擇候選角點并在這些候選角點中選擇L-曲線的最佳角點。

2.2 分步選取法

前面提出的奇異值分解求解病態方程組的正則化因子的分步選取方法，主要分為兩個階段，①在L-曲線上抽取不同數目的點構建新的L-曲線，尋找每條新L-曲線的最佳角點作為候選角點，可構成一個候選角點列表；②從第一階段得到的候選角點列表中選出最佳角點，在候選角點列表中，最佳角點應該是到達垂直分支之前的最后一個候選角點。如果候選角點分布于L-曲線的垂直分支和水平分支上，那么最佳角點為L-曲線到達垂直分支前的一個候選角點；如果曲線的垂直分支上沒有任何候選角點，那么最后一個候選點為最佳角點。

2.2.1 搜索候選角點

在算法的第一階段中,基于新的L-曲線搜索候選角點主要包括兩個準則:①基于角度搜索候選角點;②基于最小距離搜索候選角點。

1)基于角度搜索候選角點：對于L-曲線上n個離散點按照λ減小的順序依次編號為Pi(i=1,2,…n)，用向量vi,i+1(i=1,2,…,n-1)表示點Pi與點Pi+1構成的向量。首先計算向量的模長并對向量作歸一化處理使得‖vi,i+1‖2=1，再選擇其中模長較長的p(p=min(5,n-1))個向量，這些向量構成新的L-曲線，計算相鄰向量之間的叉積

delta=(vx,x+1)1(vy,y+1)2-(vx,x+1)2(vy,y+1)1

(8)

其中：x

2)基于最小距離搜索候選角點：定義水平向量v=(-1,0)T，對于上述子L-曲線上的p個向量，計算向量與水平向量v之間的叉積，公式下角標的1和2表示向量的第1個和第2個坐標

-(vj,j+1)2

(9)

|Δ|的大小反映歸一化向量在垂直方向上分量的大小，在這些向量中選擇垂直分量最小的向量vlh,lh+1作為具有水平特征的向量，選擇垂直分量最大的向量vlv,lv+1作為具有垂直特征的向量，但具有垂直特征的向量必須位于水平特征向量左側。那么“原點o”可以定義為Plh和所處的水平線與向量vlv,lv+1所在的直線的交點，然后在整個L-曲線上的所有離散點中選擇與該原點o具有最小歐幾里得距離的點作為候選角點[10]。

2.2.2 選定最佳角點

判斷第一階段得到的候選角點列表，如果列表為空，那么認為該L-曲線不存在最佳角點，否則將點P1添加到候選角點列表中，并對列表按照索引值由小到大排序，假設原來列表中候選角點個數為k，那么它們構成k個向量vj,j+1(j=1,2,…,k)，對向量作歸一化處理，判斷每個向量的垂直變化量與水平變化量,即正則化解范數的變化量與對應的殘差范數的變化量，如果正則化解范數變化量不小于相應的殘差范數變化量，那么認為該向量具有垂直特征，否則認為該向量具有水平特征。如果所有的向量都是具有水平特征，那么選擇候選角點列表中的最后一個候選角點作為L-曲線的最佳角點，如果存在具有垂直特征的向量，那么選擇到達垂直特征向量前的最后一個候選角點作為最佳角點。

3 數值實驗分析

3.1 病態問題分析

高階希爾伯特(Hilbert)矩陣方程組具有嚴重病態特征，為測試上述方法求解病態線性方程組的能力，采用100階希爾伯特方程組為例進行測算。希爾伯特矩陣方程組如式(10)所示。

Ax=b

采用Oligo 7.0軟件設計引物。針對SSR位點查找的結果，在二、三、四、五、六重復單元的位點中分別設計出10對引物，共設計出50對引物，委托上海生工生物工程股份有限公司進行合成。通過瓊脂糖凝膠電泳和變性聚丙烯酰胺凝膠電泳篩選出具多態性的引物。篩選獲得的多態引物將在其5’端添加羥基熒光素標記（FAM），進行熒光引物的合成。

(10)

圖2 100階希爾伯特矩陣的奇異值分布圖Fig.2 Singular value distribution of the 100-order Hilbert matrix

100階希爾伯特方程組經奇異值分解后的主對角元素如圖2所示，從圖2中可以看出，奇異值元素非負且降序排列，曲線首支下降較快，尾支曲線近于水平且奇異值很小，病態問題非常嚴重，如果不對奇異值進行處理，其求解效果幾乎無異于其它常規方法。通過給定不同的正則化因子可以削弱方程的病態程度，正則化因子λ通過λ(0)=1,λ(i)=λ(i-1)/100,i=1、2、…、9給出，方程組的求解誤差以最大絕對誤差ξmax=max(|1-xi|),i=1、2、…、100為評判標準。

圖3 不同正則化因子對應的最大絕對誤差曲線Fig.3 Maximum absolute error curve corresponding to different regularization factors

圖4 不同正則化因子的100階希爾伯特方程組的解Fig.4 Solution of the 100-order Hilbert equations with different regularization factors

圖5 L-曲線圖Fig.5 L-curve

圖6 誤差曲線與 L-曲線的曲率Fig.6 Error curve and curvature of the L-curve

圖7 曲線的拐點區域局部放大曲線Fig.7 Partial magnification curve of the inflection point of the L-curve

給定不同的正則化因子并逐一求解希爾伯特方程組，解的誤差表現出明顯變化，如圖3所示。從圖3中可以看出,隨著正則化因子減小，解的最大誤差逐漸減小并達到最??；當正則化因子繼續減小時，解的最大誤差又逐漸增大；當正則化因子減小到很小如1E-18時，相當于沒有做奇異值處理。給定不同的正則化因子并逐一求解希爾伯特方程組，其對應的解如圖4所示。從圖中可以看出，當正則化因子取為1E-4～1E-6時幾乎接近精確解，解的最大誤差基本在1E-2～1E-3內變化。

3.2 正則化方法的測試與分析

3.2.1 曲率法計算結果

利用文獻[6]給出的最大曲率法對圖5的L-曲線計算曲率，繪制誤差曲線與曲率曲線疊合如6所示，可以看出,平均均方誤差最小時對應的正則化因子在1E-12～1E-11之間，而最大曲率的對應的正則化因子約為1E-15，兩者不重合，說明曲率法計算的正則化因子不是“最佳”正則化因子，為此需進一步分析原因。

取圖5拐點區的數據重新繪制曲線如圖7所示，數據區的正則化因子在1E-10～1E-16范圍內變化。可以看出，在L-曲線的斜率突變區，曲線并不平滑，導致樣條曲線擬合可能存在誤差，使得計算得到的最大曲率點出現偏差，沒有對應“最佳”正則化因子。

表1 候選角點的坐標(ρ,η)和正則化因子λ

圖8 逐步尋優法搜索最佳角點的結果Fig.8 Results of searching the optimal corner point by stepwise optimization method

圖9 L-曲線拐點處放大的結果Fig.9 The result of amplification at the inflection point of the l-curve

3.2.2 分步選取法計算結果

第一步：在L-曲線上選取候選角點。L-曲線上的30個離散點構成29個向量，選取模長最長的5個向量，根據相鄰兩向量間的角度大小關系判斷是否存在候選角點，并將搜索得到的候選角點添加至列表；增加選取的向量的個數，重復判斷角度大小關系來搜索候選角點。候選角點的坐標如表1所示，候選角點的位置如圖9中紅色圓點標識所示。

第二步：在候選角點中選取最優點。表中所示5個候選角點與第1個離散點構成5個向量，對向量作歸一化處理，判斷每個向量的垂直變化量與水平變化量的大小關系，如果垂直變化量大于水平變化量，那么認為該向量具有垂直特征，否則認為該向量具有水平特征。如果所有的向量都是具有水平特征，那么選擇候選角點列表中的最后一個候選角點作為L-曲線的最佳角點，如果存在具有垂直特征的向量，那么選擇到達垂直特征向量前的候選角點作為最佳角點。從5個候選角點中搜索最佳角點為(ρ,η) = (-63.670 2, 4.605 17)，即對應著最佳正則化因子，圖8中紫色方塊為最佳正則化因子的位置。由圖8和圖9可以看出，分步選取法搜索最佳角點的效果較好。

4 結論

通過對基于奇異值分解法求解病態線性方程組的正則因子分步選取方法的研究，得到以下幾點結論：

1)奇異值分解法是求解病態線性方程組的最有效的方法之一，通過給定最佳正則化因子，可以使病態線性方程組的求解結果接近真實解。

2)基于L-曲線給出了一種給定最佳正則化因子的分步選取方法，從測試結果來看，結果較好，并且優于最大曲率法。

3)可以將分步選取方法與奇異值分解法相結合，用于求解小規模的地球物理反演問題。對于大規模地球物理反問題，由于奇異值分解法的分解過程的耗費時間與方程階數呈指數遞增，在加快計算效率方面,仍有較大的改進空間。