正則化HSS預(yù)處理鞍點(diǎn)矩陣的多尺度算法

董朝麗，汪鈺斌

(江西農(nóng)業(yè)大學(xué)南昌商學(xué)院，江西 共青城 332020)

1 引言

在眾多重要領(lǐng)域，例如計(jì)算流體力學(xué)、電磁場推算及對眾多偏微分方程的有限差分與有限元離散中，均會生成大規(guī)模的線性方程組[1]。鞍點(diǎn)問題作為一種特殊的線性方程組，在線性彈力學(xué)、圖像處理等科學(xué)研究中擁有極為廣泛的應(yīng)用[2]。

近年來，相關(guān)研究學(xué)者對種類繁多的鞍點(diǎn)求解問題進(jìn)行了不同方面的研究，并給出很多行之有效的求解方法。曹陽等[3]提出正則化HSS預(yù)處理鞍點(diǎn)矩陣的特征值估計(jì)，分析復(fù)特征值及實(shí)特征值的上下邊界，并以特征值均為實(shí)數(shù)作為充分條件，最終證明該方法測得的復(fù)特征值更精確。董貝貝等[4]提出求解鞍點(diǎn)問題的廣義正定和反Hermitian分裂方法，利用矩陣的正定分裂構(gòu)造鞍點(diǎn)矩陣的兩種分裂格式，結(jié)合迭代收斂充要條件，利用分裂格式構(gòu)造算法迭代。劉李楠等[5]提出基于超松弛預(yù)處理的精細(xì)積迭代反演算法，將方程組求解歸結(jié)為初值問題的極限形式，利用超松弛法預(yù)處理正演所得鞍點(diǎn)矩陣，降低條件數(shù)以減少測量擾動對反演計(jì)算的影響。鞍點(diǎn)系統(tǒng)中系統(tǒng)矩陣對角塊具有奇異性，在系數(shù)矩陣階數(shù)較大的情況下，受到計(jì)算機(jī)CPU、內(nèi)存與偏差等眾多外界因素干擾，上述的某些鞍點(diǎn)求解方法就不能給予較為準(zhǔn)確的收斂解。

為了完善鞍點(diǎn)問題的計(jì)算收斂速度，減少計(jì)算機(jī)CPU占用率，使計(jì)算結(jié)果更加精準(zhǔn)，本文提出一種正則化HSS(Regularrized Hermitian and skew-Hermitian splitting，RHSS)預(yù)處理鞍點(diǎn)矩陣的多尺度算法。在鞍點(diǎn)求解時(shí)，加入正則化方法，令方程的收斂效率得到顯著提高，同時(shí)對HSS方法進(jìn)行預(yù)處理，得到全新的預(yù)處理子，且對其參變量采取擇優(yōu)選取，保證算法的可靠性。通過多尺度參變量數(shù)值仿真，結(jié)果證明本文算法的收斂性極強(qiáng)，具有較高的實(shí)用價(jià)值。

2 基于正則化選取最優(yōu)參數(shù)

因?yàn)榉匠滔禂?shù)矩陣具備相當(dāng)程度的復(fù)雜度，只使用單純的迭代算法來求解式(1)，其計(jì)算效率較低。而正則化技術(shù)則可以充分改進(jìn)系數(shù)矩陣的特質(zhì)，令方程的求解過程更加便利。

下文將求解式(1)的大規(guī)模稀疏線性代數(shù)方程組

Ax=b

(1)

式(1)中，A為非奇異矩陣，x，b為列向量。在使用正則化手段的過程中，挑選適當(dāng)?shù)恼齽t化系數(shù)是重要的一步，如果選擇的系數(shù)太大，即便能夠顯著改善矩陣條件數(shù)，也會使獲得的結(jié)果與真實(shí)結(jié)果相差較多，準(zhǔn)確度不高；若系數(shù)太小又會令條件數(shù)無法得到改善，方程收斂速率也很小。因此找到一個(gè)合適的正則化系數(shù)是加快方程收斂效率的重要保證。

Tikhonov方法是使用最為普遍且高效的正則化方法[6]，其Gauss Markon模型可表示為

(2)

對式(2)采用信道估計(jì)算法可得到如下方程

(3)

若式(3)中的條件數(shù)較大，也就是在式(3)的方程為病態(tài)時(shí)，該方程求解結(jié)果準(zhǔn)確性較低，因此獲取的參數(shù)估值也相對不可信。按照Tikhonov正則化原則，可把式(2)的正則化估計(jì)標(biāo)準(zhǔn)描述為

(4)

式(4)中，K為一個(gè)對稱的正則化矩陣。由此，可將式(2)方程的求解過程重新表達(dá)為

(5)

在K=I的情況下，I為單位陣，可通過式(5)得到

(6)

與式(3)相比較，式(6)左端系數(shù)陣增添了α、I，如果α≥0，則稱其為正則化參變量。因?yàn)樵鎏砹诉@兩項(xiàng)系數(shù)可以預(yù)防方程系數(shù)陣產(chǎn)生病態(tài)問題，所以可獲得參變量的穩(wěn)定解?？梢钥闯?，在均方誤差的跡為最小的情況下，Tikhonov正則化方法可以準(zhǔn)確得到較為正確的結(jié)果，提升了算法的可靠性。由此，可將式(1)方程利用Tikhonov正則化方法轉(zhuǎn)變成求解最小平方解的問題，具體表達(dá)為

(7)

式(7)中，L為一個(gè)和系數(shù)矩陣奇數(shù)相等的矩陣，λ為正則化參變量，且此參變量是正數(shù)。一般情形下，假設(shè)L是單位矩陣I，則式(7)與(8)的求解方程為相等關(guān)系

(ATA+λ2I)xλ=ATb

(8)

式(8)中AT是系數(shù)矩陣A的轉(zhuǎn)置，xλ是使用正則化處理后的值。在求解式(8)時(shí)，如果系數(shù)λ的值過大，那么就會導(dǎo)致系數(shù)矩陣ATA+λ2I和初始方程的系數(shù)矩陣ATA相差較多，最后得到的矢量值與真實(shí)矢量值有較大誤差；相反，如果系數(shù)λ的值過小，正則化算子ATA+λ2I和初始方程ATA過于靠近，就無法具備加速收斂的效果。并且，在式(8)方程的真實(shí)操作下，ATA會致使運(yùn)算數(shù)量成倍上升，條件數(shù)的數(shù)量也會增多，令方程的收斂過程變得更加困難，所以，應(yīng)該按照原始方程系數(shù)矩陣的特征來適時(shí)調(diào)整正則化方程[7]。

根據(jù)ATA算子會引發(fā)條件數(shù)增多這一規(guī)律，在此次研究方法中使用ATA+λ2I算子并不是一個(gè)最好的選擇。此時(shí)，可把式(8)轉(zhuǎn)換成等價(jià)模式

(Α+λI)xλ=b

(9)

式(9)將式(8)方程內(nèi)的λ2利用正實(shí)數(shù)λ替換。能否選擇一個(gè)合適的正則化參變量λ是求得精確結(jié)果的重要依據(jù)。通過上述方程式可知，λ越大，就會越趨近主對角占優(yōu)，條件數(shù)越小，方程組的求解越簡單。

在式(9)方程中，A+λI算子和原始系數(shù)矩陣的算子的根本區(qū)別就是主對角元素，所以λ的選取和系數(shù)矩陣有著密切關(guān)聯(lián)。若對額外的矩陣信息為已知狀態(tài)，如特征值、條件數(shù)等，則較容易獲取正則化參變量λ。

3 HSS預(yù)處理鞍點(diǎn)矩陣的多尺度算法

在明確如何使用正則化方法選取合適的系數(shù)后，下面對鞍點(diǎn)矩陣的HSS預(yù)處理進(jìn)行進(jìn)一步分析，以此將鞍點(diǎn)問題求解的性能實(shí)現(xiàn)最優(yōu)。

HSS是一種對稱/反對稱分裂的單參數(shù)分裂迭代方法，并驗(yàn)證了在α>0時(shí)，HSS分裂迭代無條件收斂與線性方程組的唯一解。通常情況下，可將迭代法的模式表示為

x(k)=φk(x(k-1)，…，x(k-l))，k=l，l+1

(10)

式(10)中，φk為迭代算子，x(k-1)，…，x(k-l)為迭代初始值，此種迭代方法為l步迭代算法。在l=1的狀態(tài)下，則被稱之為單步迭代算法；在φk和k互不相關(guān)時(shí)，也就是φk=φ，則稱之為定常迭代法，反之稱為非定常迭代法。本文將HSS迭代方法作如下描述：

首先對系數(shù)矩陣A進(jìn)行分裂

A=H+S

(11)

(12)

通常狀況下，迭代方法的收斂速率和線性方程的系數(shù)矩陣特質(zhì)有較大關(guān)聯(lián)[8]。在系數(shù)矩陣無法滿足相關(guān)條件時(shí)，迭代算法的收斂速率相對緩慢，此時(shí)采取預(yù)處理可以很好地解決這一問題。通過預(yù)處理后的線性方程，其系數(shù)矩陣特征被重新完善，可將處理較為困難的線性方程求解問題轉(zhuǎn)變成求出其近似解問題[9]。接下來對HSS迭代方法的預(yù)處理子采取進(jìn)一步研究。

在不更改預(yù)處理矩陣結(jié)構(gòu)的基礎(chǔ)上，為了提升鞍點(diǎn)問題的求解速度，本文將HSS預(yù)處理子的衍生模式表示為

(13)

待求解的線性系統(tǒng)自身是以塊結(jié)構(gòu)化所構(gòu)成的，鑒于該結(jié)構(gòu)的良好性能，提出一新的HSS預(yù)處理子(NHSS)

(14)

這樣，就形成了全新的兩步迭代方程組，可將其描述為

(15)

考慮到有如下預(yù)處理子的存在

Mα，β=(Σ+Σ1)-1(Σ1+H)(Σ+S)

(16)

假設(shè)μ在Σ中塊結(jié)構(gòu)處起到了正則化作用，利用相關(guān)運(yùn)算可得到式(17)，β可以選擇任意值

(17)

相應(yīng)地轉(zhuǎn)變參變量范圍也能重新組合成另一種預(yù)處理子，將其表達(dá)為

(18)

因?yàn)槠渲幸粋€(gè)選擇參變量和原始矩陣內(nèi)固定的參變量是相同的，因此能夠很直接地看出兩種預(yù)處理子的差別，兩種預(yù)處理子內(nèi)的Σ與Σ1為相等關(guān)系。

本文想要構(gòu)建一種多參數(shù)化的預(yù)處理子，此種預(yù)處理子可以不受到原始矩陣內(nèi)參變量的影響，防止參變量方位調(diào)整所造成的差異化問題。所以要讓Σ滿足以下條件

(19)

式(19)中，四個(gè)參變量均和原始矩陣A內(nèi)的參變量無關(guān)且各不相同，不是一般性，可構(gòu)建出如下預(yù)處理子

M=(Σ+Σ1)-1(Σ+H)(Σ1+S)

(20)

從式(20)中，可以得到

N=M-A

=(Σ+Σ1)-1(Σ-H)(Σ1-S)

(21)

式(20)中，N為M關(guān)于A的剩余矩陣。為了更加具體地分析預(yù)處理后的鞍點(diǎn)矩陣的特征值分布狀態(tài)，需要擇優(yōu)選取預(yù)處理子內(nèi)的參變量，以此保證算法的收斂速率。

將式(20)中的M作為例子，可以看出能夠?qū)分解成

A=M-N

(22)

N=(Σ+Σ1)-1(Σ-H)(Σ1-S)

(23)

同時(shí)在式(22)中還隱含M-1A=I-M-1N，由此可以得到

λ(M-1A)=1-λ(M-1N)

(24)

(25)

(26)

根據(jù)式(26)可以看出β對N的F-范及預(yù)處理子M的預(yù)處理效果不會產(chǎn)生任何影響。由此可得出以下結(jié)論：

正則化參變量μ對運(yùn)算過程具有十分重要的作用，若Σ—H的塊取值是0，則Σ1—S的塊對N的F-范不會發(fā)生影響，即塊內(nèi)全部的參變量可以選擇任意數(shù)值；

利用最小化N的跡找出最優(yōu)參變量，且該參變量可將Σ—H中的塊取值為0；

按照式(25)和(26)中A的有效多參數(shù)，HSS預(yù)處理子可簡化為

M=(Σ+Σ1)-1(Σ1+H)(Σ+S)

(27)

通過上述過程，可以加速鞍點(diǎn)矩陣的收斂速度，使線性方程的求解更準(zhǔn)確，保證了算法的有效性與適用性。

4 仿真研究

為了驗(yàn)證研究方法的可靠性，進(jìn)行了仿真。仿真環(huán)境采用MATLAB R2012a版軟件，文獻(xiàn)[10]實(shí)驗(yàn)中所使用IFISS軟件包生成此次實(shí)驗(yàn)所用的線性系統(tǒng)。所用計(jì)算機(jī)的處理器為Inte1(R) Core(TM)2 Duo CPU T5450@1.66GHZ，Pentium雙核2.8 GHz，運(yùn)行內(nèi)存為2GB，選用64位Windows 7系統(tǒng)。將研究方法與文獻(xiàn)[3]方法(正則化HSS預(yù)處理鞍點(diǎn)矩陣的特征值估計(jì))和文獻(xiàn)[4]方法(求解鞍點(diǎn)問題的廣義正定和反Hermitian分裂方法)進(jìn)行多尺度參變量數(shù)值比較。

若原始矢量均為零矢量，右端矢量是c=He，其中e=ones(m+n，1)，迭代數(shù)量為IT，則相對殘差可描述為

(28)

式(28)中，(x(k)T，y(k)T)T為鞍點(diǎn)矩陣的第k次迭代近似解，在殘差符合RES<10-6的條件時(shí)，則算法結(jié)束。

考慮如下鞍點(diǎn)問題，其系數(shù)矩陣依次為

(29)

在此例中，將給定常量v的值設(shè)置為0.1和0.01，，針對不同的v，取相應(yīng)的p=15，20，25，同時(shí)α、β∈[0.0001，1]，數(shù)值實(shí)驗(yàn)結(jié)果如表1、表2所示。

表1 v=0.1的數(shù)值實(shí)驗(yàn)結(jié)果

表2 v=0.01的數(shù)值實(shí)驗(yàn)結(jié)果

從表1和表2中可知，在α與β的取值不同時(shí)，不管是迭代次數(shù)還是CPU，研究方法的數(shù)據(jù)實(shí)驗(yàn)結(jié)果均比其它兩種方法要好，證明研究方法的可靠性較高，CPU占用率較低。因?yàn)閷€性方程進(jìn)行了預(yù)處理，重新完善了系數(shù)矩陣特征，所以在鞍點(diǎn)矩陣求解有較好的收斂性，使數(shù)值結(jié)果的準(zhǔn)確性較高。

5 結(jié)論

為了加強(qiáng)線性方程中鞍點(diǎn)問題求解的收斂性，提高求解的準(zhǔn)確性，本文提出一種正則化HSS預(yù)處理鞍點(diǎn)矩陣多尺度算法。該方法可以顯著提高求解過程的收斂性，且CPU占用率較低，擁有較好的適用性，可為大型線性方程求解提供新的思路。